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ABSTIACT 


This  paper  is  based  on  a  paper  by  Beylkin  in  which  a  leading  order 
asymptotic  theory  for  inversion  of  acoustic  data  is  presented.  The  method 
is  based  on  earlier  work  by  Beylkin  in  the  theory  of  pseudo-differential 
operators  and  generalized  back  projections  or  Radon  transforms.  The  back 
projection  or  inversion  is  carried  out  with  respect  to  a  general  [c(z.y,z)] 
background  sound  speed.  The  asymptotic  limit  of  interest  is  high  frequency. 
The  inversion  operator  is  given  as  an  integral  of  the  observed  data  over 
frequency  and  over  the  observation  surface.  Beylkin  claims  that  his  result 
is  useful  for  finding  discontinuities  in  the  sound  speed,  but  he  does  not 
make  clear  how  this  is  to  be  done  in  practice.  I  show  how  to  modify 
Beylkin' s  inversion  operator  to  obtain  an  operator  whose  output  is  an  array 
of  singular  functions,  one  for  each  reflector  (discontinuity  surface  of  the 
sound  speed)  in  the  subsurface.  The  singular  function  of  a  surface  is  a 
Dirac  delta  function  whose  support  lies  on  that  surface.  Thus,  the  array  of 
singular  functions  produces  a  reflector  map  of  the  subsurface.  The 
validity  of  modification  of  Beylkin's  inversion  operator  is  verified  by 
applying  it  to  band  limited  Born-approximate  and  then  Eirchhof f-approximate 
representations  of  the  upward  propagating  wave  field.^  Mul ti -dimensional 
stationary  phase  is  applied  to  the  spatial  integration  oveK^the  variables  of 
the  field  representation  and  the  variables  of  the  observation  surface.  It 
is  confirmed  that  the  output  is  proportional  to  the  band  limited  singular 
functions  of  the  reflectors  and  further  that  one  can  estimate  the  jump  in 
velocity  across  each  reflector  from  the  peak  amplitude  of  the  output  on  each 
reflector.  This  is  done  for  the  cases  of  common  (or  single  fixed)  source, 
common  receiver,  and  common  (or  fixed)  offset  between  source  and  r^eiver, 
with  zero  offset  or  backscatter  as  a  special  case  of  the  last  of  th^e. 


A(x,Xs) 


a(x) 
c  (x) 
D(£,b)) 

65(s) 

F(u)) 


Amplitade  of  ray~theoretic  or  WEBJ  Green's  function  for 
background  sound  speed  with  source  at  x^  and  observation 
point  X*  See  Appendix  A. 

Perturbation  in  sound  speed.  Eq.  (3). 

Reference  sonnd  speed.  Eq.  (3). 

Observed  data,  ug(Xj.(£) ,  x^Ct)  »•*>)  •  Eq.  (5). 

gX,))  Band  limited  Dirac  delta  function  of  4  (defined  below) 
with  X',  2^  and  x^  evaluated  at  the  stationary  point,  defined 
by  Eq.  (Cl),  as  functions  of  x* 

Band  limited  Dirac  delta  function  of  s,  arc  length  normal  to 
a  surface  on  which  s  =  0.  Band  limited  singular  function  of 
the  surface.  Eq.  (38). 

Jump  in  a(x)  across  the  surface  Sj.  Eq.  (24). 

Filtered  (smoothed  and  tapered)  source  function  in  the 
Fourier  domain. 

[j.)  Phase  of  inversion  operator  applied  to  Born-approximate  or 
Kirchhof f-approximate  field  data.  Eq.  (9). 

Hessian  matrix  of  the  phase  of  the  inversion  operator  applied 
to  Born-approximate  or  Kirchhof f-approximate  field  data. 
Eq.  (33). 

First  fundamental  form  of  differential  geometry  evaluated  on 
the  reflector  S-.  Eq.  (29). 


First  fundamental  foon  of  differential  geometry  evaluated  on 
the  source /receiver  surface  evaluated  at  x^.  Eq.  (B9). 


Yj  (x) 


YjB<5) 


j 


Singular  function  of  the  surface  Sj. 

Band  limited  version  of  yjCx). 

Determinant  defining  transformation  from  variables 
to  Fourier  wave  vector  k.  Eq.  (6) . 

Jacobi  determinant  of  ray  theory.  Eq.  (A9). 


Approximate  Fourier  variable  of  the  inversion  theory. 
Eq.  (11). 


Upward  normal  vector  on  reflectors. 


s 


S. 

Sj.  J  >  1 


o 


o  =  (Oji^.O,) 


VtCx.Xj.).  Eq.  (A5). 

Kay-theoretic  reflection  coefficient.  Eq.  (47). 

Reflector  in  Kirchhoff  representation  of  upward  propagationg 
wave. 

Observation  surface. 

Reflectors.  surfaces  of  discontinuity  of  o(x)  the 

subsnrf  ace. 

Domain  of  integration  in  ^-variables. 

Running  parameter  along  rays.  Appendix  A. 

Parameters  used  to  define  a  reflecting  surface. 


Ray- theoretic  travel  time  between  x  and  x^. 


See  Appendix  A. 


Angle  between  the  noimal  to  a  surface  at  the  point  x 
ray  from  or  z^  to  z'>  under  the  stationarity  conditions 
(Cl).  Opening  angle  between  these  rays  and  nomal  subject  to 
Snell's  law  of  reflection. 

Observed  field  data.  Eq.  (5). 

*  s 

'(z)  Soimd  speed.  Eq.  (2). 

Point  at  which  the  output  of  inversion  operator  applied  to 
D(£>o))  is  to  be  evaluated. 

=  z'(a)  Point  on  reflecting  surface. 

L  ,  x.  Source  and  receiver  coordinates,  repsectively.  Eq.  (1). 


Parameters  labelling  source  and/or  receiver  points;  i.  e., 
5,  =  and/or  z^  =  z^  (5)*  Eq.  (1). 


1.  mnoDocrioN 


This  paper  is  based  on  the  brilliant  paper  by  Beylkin  [1985].  In  that 
paper.  the  anthor  presented  a  theory  for  asymptotic  inversion  of 
observations  for  the  acoustic  wave  equation  to  estimate  discontinuities  in 
the  sound  speed.  The  method  allows  for  an  assortment  of  possible 
source /receiver  configurations,  broad  enough  to  accomodate  most  of  the  cases 
of  interest  in  seismic  exploration  and  other  applications.  For  example,  the 
method  applies  to  zero-offset  data;  common  (or  single)  source,  multi¬ 
receiver  array  data  (or  the  reverse);  or  fixed  offset  data.  The  inversion 
of  the  data  is  an  integral  over  the  source /receiver  array.  Extension  to 
other  wave  equations  is  also  quite  apparent  from  the  presentation. 

Beylkin' s  results  are  couched  in  the  language  of  pseudo-differential 
operators  and  generalized  back  projections  or  generalized  inverse  Radon 
transforms.  He  also  does  not  make  precise  the  manner  in  which  his  method 
actually  produces  the  discontinuities  in  the  sound  speed.  The  theory  only 
predicts  that  his  ir.tegral  solution  is  a  leading  order  high  frequency 
asymptotic  inversion  operator  without  relating  the  output  to  the  information 
being  sought. 

The  purpose  of  this  paper  is  twofold.  First.  I  make  more  precise  the 
manner  in  which  Beylkin' s  method  provides  an  asymptotic  solution  of  the 
inverse  problem.  1  start  from  one  of  Beylkin' s  interpretations  of  his 
inversion  scheme  to  view  the  output  as  a  high  frequency  Fourier  inversion  of 
band  limited  data  for  the  perturbation  in  sound  speed.  1  show  that  by 
modifying  the  integrand  of  the  operator  by  a  scale  factor,  it  will  produce 


asymptotically  an  array  of  scalod  aingmlar  fmnctlons  of  the  surfaces  of 
discontinuity  of  the  sound  speed.  These  surfaces  are  the  reflectors  in  the 
subsurface.  The  singular  function  of  a  surface  is  a  Dirac  delta  function 
whose  support  is  On  the  surface.  Thns>  knowledge  of  the  singular  functions 
is  equivalent  to  mathematical  imaging  of  the  reflecting  surfaces.  The 
scaling  factor  for  each  reflector  is  a  known  function  of  the  jump  in  sound 
speed  across  the  reflector.  That  known  function  takes  a  different  form 
depending  on  whether  one  has  used  the  Born  approximation  or  the  Kirchhoff 
approximation  to  represent  the  input  data.  In  either  case,  I  provide  a 
means  of  estimating  from  the  output  the  change  in  sound  speed  across  the 
reflector.  The  estimate  is  consistent  with  the  fom  of  the  input  data. 

Beylkin  derives  his  inversion  operator  for  upward  traveling  waves 
represented  by  their  Born  approximation.  The  "backwards  projection"  of  this 
approach  is  with  respect  to  an  assumed  known  reference  sound  speed.  I  begin 
from  such  a  representation,  as  well.  From  this  starting  point.  1  can  only 
have  confidence  in  the  accuracy  of  our  interpretation  of  the  output  of  our 
method  for  small  perturbations  in  sound  speed. 

The  output  that  is  produced  depends  on  a  certain  open  angle  between  two 
particular  rays  in  the  subsurface.  This  angle  is  not  known  because  it 
varies  from  point  to  point  in  the  subsurface.  At  each  point,  it  depends  on 
the  directions  of  incidence  and  reflection  of  a  pair  of  specular  rays  from 
an  (unknown)  source/receiver  pair  on  the  datum  surface.  I  introduce  an 
alternative  inversion  integral  which,  in  the  Born  limit  produces  an  output 
whose  peak  value  on  a  reflector  is  independent  of  that  opening  angle. 


1  then  consider  the  application  of  my  inversion  operator  to  Kirchhoff 


approximate  data  for  a  single  reflector.  For  such  a  representation.  I  need 
not  assume  small  perturbations,  but  only  high  frequency,  which  is  already 
assumed  in  this  theory.  Here.  I  am  able  to  show  that  the  weighting  on  the 
singular  function  of  the  surface  is.  to  leading  order  asymptotically,  the 
full  angular  dependent  reflection  coefficient,  independent  of  the  size  of 
the  jump  in  velocity  across  the  reflector.  Thus,  it  would  seem  that  what 
needs  be  "small"  in  this  formalism  is  the  error  between  the  background  sound 
speed  above  a  given  reflector  and  the  "true”  sound  speed  above  the 
reflector. 

The  angle  of  the  reflection  coefficient  is  as  described  above.  I  show 
how  this  angle  can  be  determined  by  exploiting  the  two  inversion  integrals 
already  proposed.  Thus,  in  terms  of  numerical  processing  one  need  compute 
only  one  additional  sum  with  summand  given  in  terms  of  previously  computed 
elements.  From  these  two  outputs,  one  obtains  a  reflector  map.  an  angularly 
dependent  reflection  coefficient  and  an  estimate  of  the  (cosine  of  that) 
angle. 

I  then  propose  a  third  inversion  integral.  This  one  has  the  property 
that  its  peak  amplitude  on  the  reflecting  surface  is  equal  to  the  product  of 
the  angularly  dependent  reflection  coefficient  and  the  area  under  the 
temporal  filter  of  the  original  time  signal.  This  last  result  is  the  most 
esthetically  appealing,  but  offers  little  practical  advantage  over  the  other 
two  operators.  In  any  case,  at  least  two  inversion  integrals  must  be 
computed  to  determine  the  unknown  angle  at  each  output  point  and  then  the 
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jump  in  velocity  from  the  angnlar-de pendent  reflection  coefficient. 

This  verification  also  suggests  a  recursive  application  of  the 
inversion  formalism.  That  is,  starting  from  the  upper  surface,  each  time  a 
'major*  reflector  is  imaged,  the  background  sound  speed  is  updated  to 
account  for  the  new  information  and  data  is  processed  deeper  into  the 

section  until  a  new  major  reflector  is  imaged.  The  method  is 

pointwise ,  hence  lending  itself  to  this  type  of  recursive  impl mnentation. 

In  these  results,  the  upper  surface  is  allowed  to  be  curved.  Thus,  the 
inversion  one  produces  eliminates  two  preprocessing  steps  usually  applied  to 
seismic  data.  The  first  is  a  static  correction  for  variable  height  of  the 
sour ce/re ce iver  array.  The  second  is  stacking  to  produce  an  'equivalent* 
zero  offset  (backscatter )  data  set. 

Central  to  this  derivation  of  these  results  is  the  method  of  multi¬ 
dimensional  stationary  phase.  This  brings  me  to  the  second  point  of  this 
paper,  namely,  to  provide  a  more  classical  verification  of  the  asymptotic 
validity  of  our  modifications  of  Beylkin's  inversion  operator.  Also,  the 
interpretation  in  terms  of  imaging  of  reflectors  and  estimating  reflection 
coefficients  arises  in  a  natural  way  in  this  method.  Furthermore,  the 

method  predicts  that  such  imaging  will  occur  only  at  those  points  on  the 
reflector  for  which  there  are  a  specular  pair  of  rays  from  the  source  and 
the  receiver  to  the  surface  point.  This  ties  the  inversion  back  to  the 

required  sour  ce  /  re  ce  iver  array  necessary  for  imaging  the  reflectors  in  the 
subsurf  ace . 

I 
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The  verification  of  the  t«a-aiid-oiie~half  diaensional  (2.SD) 
Specialization  of  Beylkin's  result  has  already  been  presented  in  Bleistein. 
Cohen  and  Hagin  [198Sb].  In  that  case,  it  is  assnmed  that  the  data  is 
gathered  only  on  a  single  line  on  the  surface  and  that  all  parameters  of  the 
subsurface  are  functions  of  the  transverse  variable  along  that  line  and 
depth.  Significant  simplifications  occur  in  that  analysis  because  fourfold 
integrals  over  tvo  surfaces  appearing  in  the  analysis  here  reduce  to  twofold 
integrals  over  two  curves  in  the  2 . 5D  case.  Furthermore,  a  certain  3X3 
matrix  central  to  Beylkin's  approach,  appearing  in  the  present  analysis, 
reduces  to  a  2><2  matrix  in  the  2 . 5D  case  and  can  be  analyzed  much  more 
readily  than  here. 

The  modification  of  Beylkin  which  I  use  is  a  generalization  of  the 
method  previously  employed  in  Bleistein  and  Cohen  [1979,  1982],  Bleistein 
and  Gray  [1985],  Cohen,  Bleistein  and  Hagin  [1985],  and  Bleistein,  Cohen  and 
Hagin  [1985a]  and  [1985b].  The  essential  feature  of  this  modification  is  a 
fundamental  result  in  Fourier  analysis,  namely,  given  the  Fourier  transform 
of  a  function  with  surfaces  (3D)  or  curves  (2D)  of  discontinuity, 
multiplication  of  the  Fourier  data  by  ±ik  before  inversion,  produces  the 
array  of  singular  functions  of  the  discontinuity  curves  or  surfaces.  Here, 
k  is  the  magnitude  of  the  spatial  transform  vector  variable  and  ±1  =  sgnia 
(time  transform  variable)  in  the  present  application.  In  the  application  to 
Beylkin's  result,  I  need  only  identify  his  inversion  representation  as  an 
inverse  Fourier  transform  and  then  identify  k  f or  this  representation.  It 
then  remains  to  carry  out  the  verification  of  the  validity  of  the 
application  of  this  theory  to  upward  scattered  data  in  this  general  context. 


This  paper  contains  an  extensive  appendix  on  the  relation  between 
Beylkin's  transformation  determinant.  >  and  the  Jacobi  determinant 

which  naturally  arises  in  ray  theory.  This  analysis  is  important  if  the 
method  is  to  be  implemented  numerically.  Furthermore,  the  length  of  this 
appendix  reflects  ny  personal  interest  in  this  interplay  between  ray  theory 
and  asymptotic  inversion. 


2.  MODIFTING  BBILUN'S  KESOLT  TO  OBTAIN  TBB  SINGDLAB  FONCTION 


I  present  here  Beylkin' s  solution  to  the  aconstic  inverse  problem  in  a 
constant  density  medium.  Let  us  consider  a  seismic  experiment  carried  out 
on  the  surface  of  the  earth  in  which  the  source /receiver  pairs,  x^  and  x^, 
respectively,  are  identified  by  a  parameter  £  =  follows: 


?s  =  ?r  =  ?r^^^ 


For  example,  for  the  case  of  common  source.  Xj  would  be  a  constant  vector 
denoting  that  fixed  position  and  the  function  x^(^)  would  be  a  parametric 
representation  of  the  receiver  surface;  for  the  common  receiver  case,  the 
roles  of  x^  and  x^  would  be  reversed;  for  the  coaunon  offset  case  or  common 
midpoint  case,  both  x^  and  x^  would  vary  with 


It  is  assumed  that  uCx.Xj.oi)  is  the  response  to  an  impulsive  point 
source  at  x^  satisfying  the  wave  equation 


•  V  •  4  < 

■:-v 


•  -■ .  j 
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V’  .  •  .  ;  4  •  4 


:m-  .j 


V  n  -  ^  u  =  6  (x  -  X  )  • 

a  -  -5 


Born-approximate  inversion  is  based  on  the  assumption  that  the  propagation 
speed  can  be  written  in  terms  of  some  reference  speed.  c(x).  *nd  a  small 
perturbation  a(x)  as  follows: 


7  =  7l‘*“] 


S  •. 


The  total  field,  u  in  (2),  is  then  written  in  terms  of  an  incident  field 


which  satisfies  (2)  with  v  replaced  by  c,  and  a  scattered  field,  ,z^,<<)) , 

which  is  everything  else. 

It  is  further  assumed  that  observations  of  this  latter  field  are  made 
at  the  points  x  =  Xj.(S)  for  sources  at  the  points  x  =  Xj(£).  Beylkin' s 
inversion  formula  for  a(x)  is 


a(x) 


c*  (x) 

8jt* 


4  |i»(s»IL)| 

d  5  - 

A(x,x^)A(x,x^) 

S 

5 


dm  F(u))  exp{-i«)[T(x,x„)  +  t(x,x.)])  D(£,w)  . 

—  — s  “*  “T 


(4) 


In  this  equation  I  have  used  the  following  notation.  The  domain  of 
integration  is  the  set  of  {i^-values  which  are  required  to  cover  the 
sonrce/receiver  array.  The  notation  So  is  reserved  for  the  surface  on  which 
x^  and  x^  are  located.  The  domain  of  integration  in  u  is  limited  by  the 
'filter*  F((i)).  I  take  this  function  to  be  synunetric  and  smoothly  tapering 
to  zero  at  the  ends  of  its  support.  I  think  of  FCw)  as  a  smoothly  tapered 
version  of  the  source  wavelet.  In  (2),  there  should  have  been  some 
frequency  domain  source  function  on  the  right  side  but  I  omitted  it  there  in 
order  not  to  introduce  two  functions,  one  for  the  source  and  another  for  the 
smoothing  filter.  Both  are  contained  in  this  one  function  in  the  inversion 
formula.  The  functions  T(x,Xg)  and  A(x,x^)  It(x,Xj)  and  A(x,Xj)]  are  the 
WKBJ  or  ray-theoretic  phase  and  amplitude  of  the  Green's  function  with 


source  at  and  observation  point  z.  discussed  in  Appendix  A.  The 

function,  D(£t(<>)  is  a  shortened  notation  for  the  observed  data  at  the  upper 
surf  ace  : 

D(§,(i))  =  u„(x  (^)  ,x  (£)  ,<ii)  .  (5) 

O  “i  ~S 

The  function.  h(x.^)>  is  the  essential  element  of  Beylkin's  result.  It  is 
the  determinant. 


V[t(x,x  )  +  t(x,x  )1 
”  ”s  r 

h(x.4)  =  det  •|^V[t(x.x  )  +  tCx.x  )]  (6) 

$  r 

^(x.x  )1 

-  -s  -  -r 

It  is  assumed  throughout  that  h  ^  0  and  is  finite.  For  four 
configurations  of  and  x^C^)  of  interest  in  geophysical  experiments,  we 

show  in  Appendix  B  that  this  is  equivalent  to  the  assumption  that  there  are 
no  caustics  in  the  ray  families  between  the  output  point  x  and  the  upper 
surface.  , 

I  will  digress  here  to  provide  an  interpretation  of  the  result  (4). 
This  interpretation  is  a  synthesis  of  Beylkin’s  own  presentation  and  the 
discussion  to  be  found  in  Cohen,  Bleistein  and  Bagin  [1985].  Let  us 
consider  the  high  frequency  Born  approximation  of  the  data.  in  terms 

of  the  perturbation,  a(x) .  That  data  is  given  by 


D(§,u)  =-(!»*  c  *(x')  a(x')  A(x',x^)  A(x',x^) 


»xp{iw[T(x' ,x^)  +  t(x',x^)]}  d  x' 


This  result  is  inserted  into  (4)  to  produce  the  following. 


a(x)  ~  - 


:*U)  ,  |h(x.£)| 

- -  d  ?  - 

8n  A(x.x^)A(x,Xj.) 


dm  F((i)) 


c  *(x')  a(x')  A(x’,x  )  A(x',x  )  d*x' 


sxp  {iw  l(x,x',x^,Xj.))  • 


In  this  equation. 


f(x,x',x  ,x)  =  x(x',x  )  +  t(x',x,)  -  [t(x,x^)  +  t{x,x^)] 

-  -  -s  ~t  -  “s  -  ~r  ~  “S  “  r 


is  the  difference  of  travel  times,  source  to  input  point  to  receiver  minus 
source  to  output  point  to  receiver. 


Beylkin's  theoretical  approach  to  the  asymptotic  inversion  predicts 
that  the  dominant  critical  point  of  the  integral  (8)  is  the  point  where  i 
vanishes,  namely,  where  x'  =  x,  for  any  choice  of  Furthermore,  in  a 

neighborhood  of  that  point,  linear  approximation  (or  first  term  of  the 


- 

r  w 


a  %.  ■- 


•‘.v’.V  K 

r»r-  -1 


■■  •  .*  •*  *  -  w  _  W  a  _  •  ^  •  «  •  — _  "f  ■_.*4  I 
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Taylor  series)  for  i  yields 


f (x.x' a  VIt(x,x^) 


+  t(x,x^)] 


(x*  -  x) 


(10) 


With  this  approximation,  the  phase  in  (8)  has  the  form  of  a  Fonrier  phase, 
ik*(x'  -  x) ,  with  the  wave  vector  k  given  by 


k  =  (0  VtT(x,x^)  +  t(x,x^)]  . 


(11) 


I  view  (11)  as  defining  a  change  of  variables  of  integration  from  (u,^) 
—  three  variables  —  to  k  —  also  three  variables.  The  function,  w*h(x,^) 
is  Jnst  the  Jacobian  of  that  transformation,  so  that 

(o*|h{x,t)  I  du)d*?  =  dk*.  (12) 

Now  the  integral  in  (8)  is  seen  to  be  of  the  form  of  a  forward  and  inverse 
Fourier  transform,  producing  the  integrand  evaluated  at  x'  =  x.  That 
evaluation,  indeed,  yields  a(x)* 

This  ends  the  digression. 


It  is  this  interpretation  as  a  Fourier  integral  which  will  allow  me  to 
deduce  a  representation  for  the  reflectivity  function  from  the  solution  (8) 
for  o(x).  I  base  that  result  on  a  theory  for  identifying  surfaces  of 
discontinuity  of  a  function  from  large  jk|  band  limited  data  for  the 
function  [Cohen  and  Bleistein,  1979,  Bleistein,  1984].  The  result  is  that 


dn 


< - >  K'^ignu)  |k|  3(k)  . 


(13) 


In  this  equation,  da/dn  denotes  the  upward  normal  derivative  of  a  at  the 
surfaces  of  discontinuity  and  the  correspondence  means  that  if  we  have  the 
Fourier  data  3(k) >  obtain  the  Fourier  data  for  da/dn  by  the  indicated 
multiplication.  In  fact,  I  have  now  interpreted  (8)  as  a  Fourier  integral 
with  k  defined  by  (11).  Thus,  to  obtain  da/dn,  I  need  only  multiply  the 
integrand  in  (8)  by  the  factor 

i(signu)|k|  =  iu)  |VtT(x.r^)  +  T(x,x^)lj  .  (14) 


Let  ns  suppose  that  o(x)  has  a  discontinuity  surface,  S,  Then  the 
upward  nonnal  derivative  at  any  point  on  S  is  proportional  to  a  Dirac  delta 
function  of  arclength  along  a  curve  normal  to  S  with  peak  of  that  delta 
function  on  S  itself.  This  delta  function  is  the  singular  function  of  the 
surface  S,  It  is  through  this  depiction  of  the  singular  functions  or  the 
discontinuity  surfaces  of  a  that  processing  for  da/dn  provides  a  reflector 
map.  Furthermore,  for  each  surface  the  multiplier  of  the  singular  function 
is  equal  to  the  upward  jump  in  a(x)  across  the  surface,  namely. 


a(x-)  -  a(x+)  s  c  (x) 


v*(x-)  v*(x+) 


r(x  +  )  -  v(x~) 
2c  (x) 


(15) 


Here,  a(x~)  is  the  limiting  value  of  a  from  points  above  the  surface  and 
a(x+)  is  the  limiting  value  of  a  from  points  below  the  surface.  I  have 
written  approximate  equalities  in  these  equations  and  retained  only  the 


-  12  - 


inversion  theory  at  this  point.  That  will  be  remedied  in  Section  5  below; 
we  fflnst  content  onrselves  with  this  result  for  the  present.  One  can  see, 
however,  that  the  singular  function  theory  combined  with  the  Born,  high 
frequency  inversion  theory  provides  a  reflector  map  and  a  basis  for 
parameter  estimation. 


The  modification  (13,14)  is  now  used  in  the  expression  (4)  for  a  to 


‘  *-•  * 

>  V 


*-  V 


3 .  ASYMPTOTIC  ANALYSIS  -  FBELIMINABY  BESOLTS 


For  the  inversion  operator  {16),  I  use  for  the  Born-approximate 

data  (7)  to  obtain  an  expression  analogous  to  (8),  expressing  the  output, 
9a/dn,  in  terms  of  the  perturbation,  a,  itself.  The  result  is 


jisi^ 

-  ,  .^V  .  s  . 


■'-*  V 


V  V 


da  c 


*(x)  |h(x,^)|  |Vt(x.x  )  +Vt(x.x  )| 

__  d  4  - ^ - 

8n  A(x,x  )A(x,x  ) 

J  Jo  s  r 


id)  du  F((i)) 


(x*)  a(x')  A(x',x.)  A(x',x_)  exp  (iw  #(x,x' ,x  ,x  ) )d  x'. 

S  X  o  * 


The  phase,  #,  is  given  by  (9);  A(x,x^)  expCiwx (x, x^) }  is  the  ray-theoretic 
or  WKBJ  Green's  function  for  source  point  at  xs  and  observation  point,  x, 
with  a  similar  description  for  the  other  amplitudes  and  traveltimes;  the 
Jacobian,  h(x,^)  is  defined  by  (6). 


:m‘ 


objective  is  to  show  that  the  multi-fold  integral  in  (17) 
asymptotically  produces  an  array  of  scaled  singular  functions  of  the 

surfaces  of  discontinuity  of  a(x) .  To  do  this,  I  will  apply  the  method  of 

stationary  phase  [Bleistein,  1984]  to  the  five-fold  integral  in  x'  and  £ 

under  the  assumption  of  "high  frequency",  and  then  analyze  the  remaining 

u)-integral.  I  remind  the  reader  that  by  high  frequency  1  mean  that 


=  »  1. 
c  c 

0  0 


In  this  equation,  L  is  a  'typical  length  scale”,  c,  is  a  local  estimate  of 


the  background  sound  speed  and  f  is  the  frequency  in  Hz.  The  dimensionless 
parameter.  A,  then,  measures  typical  lengths  in  units  of  the  wave  length. 
In  practice,  a  value  of  A  of  3  is  sufficient  for  high  frequency  asymptotic 
analysis  to  adequately  approximate  a  Fourier  type  integral,  such  as  the  one 
we  have  here. 

The  integral  (17)  does  not  even  converge  unless  we  make  some  assumption 
about  the  nature  of  a(x) .  The  reason  is  that  the  amplitudes  A^z',z^)  and 
A(x',Xj)  each  decay  only  as  l/|x*|  at  infinity,  while  d*x'  is  0 (  |x '  |  *d  |x '  | ) 
at  infinity.  Thus,  to  avoid  convergence  problems  at  infinity,  I  will 
assume  that  a{x')  has  finite  support;  that  is,  this  function  is  zero  outside 
of  some  finite  domain.  This  is  not  a  serious  constraint.  Real  data  is 

always  of  finite  extent  both  spatially  and  temporally.  In  a  model 

experiment  over  all  space  and  time,  one  simply  models  this  finiteness  of 

data  as  being  a  consequence  of  the  finite  support  of  the  function  a(x').  Of 

course,  our  output  may  contain  "artifacts"  which  arise  from  the  boundaries 
of  the  input  data  set.  We  must  be  alert  to  these  and  reject  them. 

Alternatively,  one  could  taper  the  data  set  both  spatially  and 
temporally  to  minimize  the  effects  of  the  abrupt  boundaries.  Equivalently, 
I  assume  that  a(x')  vanishes  smoothly  at  the  boundaries  of  its  support 
domain,  while  I  will  still  allow  a(x’)  to  have  discontinuities  inside  that 
domain.  In  this  manner,  the  mathematical  asusmption  that  a(x')  has  finite 
support  and  vanishes  smoothly  at  the  boundary  of  its  support  is  seen  as  a 
physically  reasonable  assumption,  while  providing  a  useful  mathematical 
constraint  for  our  further  analysis. 
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Let  us  now  consider  applying  the  method  of  stationary  phase  to  (17)  in 
the  variables  x'  and  Then  the  first  derivatives  of  i  with  respect  to  all 
of  these  variables  would  have  to  be  equal  to  zero.  I  claim  that,  except,  in 
unusual  (pathological)  cases  this  cannot  happen.  More  specifically,  the 
derivatives  with  respect  to  all  of  the  x'  variables  cannot  all 
simultaneously  be  zero.  To  see  why  this  is  so,  I  write  down  those 
derivatives,  using  (9): 

V’Kx.x'.x^.x^)  =  V't(x'.x^)  +  V'r(x'.x^)-  (19) 

In  this  equation,  V'  denotes  the  gradient  with  respect  to  x’.  The  reader  is 
reminded  that  V't(x',X5)  [V't(x',x^)]  is  the  tangent  vector  to  the 
geometrical  optics  ray  from  x^  f?r^  background  sound  speed  c(x')  . 

In  order  for  Vi  to  be  zero,  then,  the  ray  directions  would  have  to  be  anti- 
col  inear.  Equivalently,  both  rays  would  have  to  be  segments  of  a  single  ray 
which  emanates  from  one  of  the  points  Xg,  x^*  travels  into  the  subsurface 
and  turns  back  to  the  surface  to  emerge  at  the  other  of  these  two  points. 
We  will  assume  that  Xg,  x^  and  x'  are  not  a  triple  of  points  for  which  this 
can  occur.  ON  the  other  hand,  however,  for  vertical  seismic  profiling  or 
other  tomographic-like  inversion  problems.  this  would  be  exactly  the 
stationary  point  of  interest. 

The  method  of  stationary  phase  dictates  that  when  there  are  no  interior 
stationary  points,  one  applies  integration  by  parts  (the  divergence  theorem 
in  more  than  one  dimension)  to  replace  the  integral  over  the  interior  by  an 
integral  over  certain  boundary  surfaces.  There  are  two  such  types  of 
boundaries:  the  actual  boundary  of  the  domain  of  integration  (where  we  have 
assumed  that  a(x')  =  0)  and  the  discontinuity  surfaces  of  the  Integrand. 


For  the  present,  I  assume  that  the  only  discontinuities  are  in  a(x')« 


itself.  At  the  end  of  Section  4.  I  will  discuss  the  effects  of 
discontinuities  in  the  background  speed  c,  as  well. 

Let  us  consider,  for  a  moment,  only  the  x'-integral  in  (17): 


I  = 


c  (x')  a(x')  A(x' 


,x  )  A(x’,x  )  exp  liu)  t(x,x',x  ,x 

~s  -  ~T  ~  “  “s  r 


)) 


d’x' 


(20) 


I  denote  by  ,  S^...,  the  discontinuity  surfaces  of  a(x').  For  simplicity, 
we  sahould  think  of  them  as  extending  across  the  support  domain  of  a(x')  and 
not  intersecting.  (This  will  simplify  our  analysis  somewhat.  It  will  be 
clear  below  that  the  case  of  intersecting  surfaces  —  lenses,  for  example  — 
can  easily  be  included  in  the  discussion.)  I  assume,  therefore,  that  a(x') 
is  given  by  different  functions  between  the  surfaces: 


'),  x' 

be  tween 

So 

and 

a(x')  = 

ai(x 

'),  x' 

between 

Sx 

and 

s,  . 

(21) 

'),  x' 

be  tween 

S 

n 

and 

^n+i' 

with  the  observation  surface  on  which  x^  and  x^  reside. 

Let  us  consider  the  integral  in  x'  over  the  domain  between  Sj  and 
where  o(x')  =  Oj(x').  In  order  to  apply  the  divergence  theorem,  we  must 
first  write  the  integrand  in  the  z'  integral  as  a  divergence.  To  do  so. 
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consider,  first,  the  identity 

Vi 


Wexp  {i(u4)  =  V 


ia)|V'#|* 


Wexp  {i(j)t) 


V' ' 


Vi 

-  W 

iujV'f  I* 


e  xp  { i  (i)i }  .  (22) 


In  this  expression,  I  have  used  W  to  denote  the  amplitude  of  the  x’-integral 
in  (20)  .  Note  that  the  second  term  on  the  right  here  is  of  the  same  form  as 
the  left  side  of  this  equation  with  a  different  amplitude  function  and  a 
factor  (id))  Thus,  asymptotically,  one  should  expect  that  the  integral  of 

this  term  will  be  of  lower  order  in  oi  [0(a)~*')]  than  the  integral  of  the  left 
side.  Therefore,  to  leading  order  asymptotically,  I  need  only  integrate  the 
first  term  on  the  right.  It  is  to  this  integral  that  the  divergence 
theoremcan  be  applied,  to  obtain 


..v.v'.v.tI 


I  = 


(x')  a.(x')  A(x',x  )  A(x',x  )  exp  {iw  f(x,x',x  ,x  )) 
j-  --S  -  -r  -  -  -s  ~T 


S.  - 

J  J+1 


(23) 


11-[V't(x'.x  )  +  V't(x'.x  )] 

_  _s  _  _j. 


V't(x',x  )  +  V't(x',x  ) 

-  “S  ~  ~T 


dS' 


In  this  equation,  ii  denotes  the  upward  normal  on  the  surface,  Sj  or  Sj^2 


This  is 

the 

outward 

normal  tc 

the  domain  of 

integration  at 

th  e  uppe  r 

surface. 

but 

it  is 

the  inward 

normal  at  the 

lower  surface. 

Thus,  one 

obtains  surface  integrals  of  opposite  sign  when  the  result  of  the  divergence 
theorem  is  expressed  in  terms  of  this  upward  normal.  In  particular,  I 
assume  that  a  =  0  on  so  that  the  integral  over  S,  is  equal  to  zero.  I 
have  also  explicitly  written  V’4  in  terms  of  the  gradients  of  the  separate 


travel  times,  from  (9). 


! 

It  is  now  necessary  to  sum  integrals  of  the  form  (23)  over  all  of  the 
separate  domains  of  definition  of  a(x').  For  each  domain  integral  we  obtain 
a  pair  of  surface  integrals.  When  this  sum  is  re-ordered  by  surfaces 
instead  of  domains,  we  obtain  a  pair  of  integrals  over  each  surface  in  which 
the  only  difference  in  the  integrand  arises  from  the  discontinuities  of 
a(x') .  Let  us  introduce  the  notation 

A  .  =  [a(x'-)  -  a(x'+)]  ,  x'  on  S.  .  (24) 

J  -  -  -  J 

As  in  the  previous  section,  (-)  denotes  the  limit  through  lower  values  of  z 
(from  above)  and  ( +)  denotes  the  limit  from  below.  Now, 


da 

da  . 

-  =5 

_ 1 

dn 

dn 

where 

(25) 


da.  c  (x) 

_ J _ 1_ 


dn  8n 


|h(x,^)|  |Vx(x,x  )  +Vt(x,x  )j 

d  5  - 


A(x,x  )A(x,x  ) 

-  -s  -  -r 


u  du)  F(u)) 


(26) 


*(x')  A.  A(x',x  )  A(x',x  )  exp  (iw  i(»,x',x  ,x 

-j--s-r  SI 


)} 


S. 

J 


n*IV'T(x',x  )  +  V't(x',x  )) 

-  s  _ -  ~r 

V’t(x',x  )  +  V't(x'.i  ) I* 

-  -S  I 


dS'  . 


In  this  integral,  the  surface  S,  is  parametrized  by  two  parameters 


It  might  be  more  proper  to  index  g  by  j,  as  well.  However,  I  will  continue 
the  discussion  below  focused  on  only  one  surface  and  hence  dispense  with 
indexing  on  this  variable.  In  terms  of  these  parameters. 


dS'  = 


wi  th  g 


Here 


(28) 


the  first  fundamental  form  of  differential  geometry  for  Sj > 


dx' 

2 

dx'  dx'  1 

-  ■  - 

X  ^ — 

det  —  *  -  1 

do^ 

do. 

do.  do  1 

1 

2 

L  k  ffl  J 

(29) 


denotes  the  vector  cross  product  and  •  denotes  the  vector  dot 


product. 


4.  ASIHPTOTIC  ANALYSIS  W  OOTFOT  FOB  BOBN-APnOXIMATB  DATA 


I  will  now  apply  the  method  of  stationary  phase  to  (26)  in  the  four 
variables  The  phase  i  is  a  function  of  these  variables  through  the 

dependence  of  x'  on  o  and  the  dependence  of  and  x^  on  Equation  (9)  is 
used  to  write  the  four  first  derivatives  of  i  in  terms  of  the  derivatives  of 
the  travel  times: 


=  V  f 

si 


t(x’ t?  )  -  t(x,x  ) 


]  •  dT  ^  ~  ] 


.  Ji^ 


*  w  •  -  ' 


*“« 

V  VV  ' 
--  ' 


[x(x-.x^)  +  x(x'.x  )  ]  ■  .  «  =  1.2. 

ffl 


The  stationary  points  in  (g.®)  s*"®  determined  by  requiring  that  these  first 
derivatives  all  be  equal  to  zero. 


In  Appendix  C>  1  discuss  the  conditions  under  which  f  is  stationary. 
The  stationary  phase  conditions  are  stated  as  equation  (Cl).  Also,  in  that 
appendix  I  show  that,  for  x  on  the  surface  Sj  for  some  fixed  j,  there  is  a 
unique  stationary  triple,  x’,  x^  and  x^,  with  x'  =  x.  This  is  shown  for  the 
following  source /receiver  configurations  of  practical  interest;  common 
source,  common  receiver  and  common  offset.  Although  I  have  only  considered 
here  the  fully  three  dimensional  problem,  this  analysis  specializes  to  the 
cases  of  2.SD  inversion. 


Sv.'A'Cv 


I  will  proceed  below  by  focusing  attention  on  this  stationary  point  on 
S.  when  x  is  in  the  neighborhood  of  S-.  That  is,  this  is  the  stationary 

J  J 

point  which  has  limit  x'  =  x  as  x  approaches  Sj.  If  there  were  no 
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source/receiver  pair  in  the  seismic  survey  under  consideration  which 
included  the  particular  and  x^  needed  to  complete  the  stationary  triple, 
then  the  asymptotic  contribution  for  that  point  x  would  be  of  lower  order  in 
w  and  almost  always  of  smaller  value  numerically  after  the  u  integration 
than  the  result  we  will  obtain  below.  Thus,  I  proceed  under  the  assumption 
that  such  a  stationary  triple  has  been  determined  and  that  the  corresponding 
values  of  a  and  ^  are  interior  points  of  the  respective  domains  of 
integration. 


StationaryPhaseEyal nation 


The  integral  in  (26)  is  evaluated  by  the  method  of  mul  ti -dimensional 
stationary  phase  in  the  four  variables  g  and  The  result  is 


i  , 

c  (x) 


c  (x') 


A(x',Xj)  A(x',x^) 

A(x,x  )  A(x,x  ) 

-  ”s  ~  -r 


|h(x,4)|  |Vt(x,|^)  +Vt(x,|j|.)| 


n-tV'r(x',x  )  +  V't(x',x  )] 
_ ~  ~s _ -  -r 

V't(x',x  )  +  V't(x',x  ) 

-  s  ~  ~r 


(31) 


In  this  equation,  gj  is  defined  by  (28)  and  Kx)  denotes  the  integral 

F(<»)  exp  {iwf  (a.x' ,x^,x^)  +  i(sign  w)  ( jt/4)  s  ig(  $^^) }  dw  ,  (32) 

This  integral,  as  well  as  the  entire  right  side  of  (31),  is  a  function  of  x, 
alone,  because  x',  x^  and  x^.  are  determined  as  functions  of  x  from  the 


stationarity  conditions,  (Cl).  The  matrix,  I§t_J ,  is  a  4X4  matrix. 


K]  = 


95,35  35.30 

i  HI  i  n 


35.30  da. da 

1  m  1  m 


t  i»iii  ^ 


det[f^f^]  denotes  the  determinant  of  this  matrix  and  sig(i^^)  denotes  the 
signature  of  the  matrix,  which  is  the  number  of  positive  eigenvalues  minus 
the  number  of  negative  eigenvalues  of  the  matrix. 


I  • 


-  ^  -  *  -  k*. 


Since  it  is  expected  that  3aj/dn  peaks  for  x  on  Sj,  «e  are  interested 
in  evaluating  the  result  (31,  32,  33)  for  x  near  Sj.  Let  us  first  consider 
the  behavior  of  the  matrix  in  (33)  when  x  is  Jn  Sj.  In  this  case,  o 

can  be  fixed  before  evaluating  the  second  derivatives  with  respect  to  5i» 
5g,.  In  that  limit,  i  ■  0;  the  entire  2X2  matrix  in  the  upper  left  hand 
corner  of  is  a  matrix  of  zeroes;  the  determinant  of  is  just  the 

square  of  the  determinant  of  the  2X2  matrix  in  the  upper  left  hand  corner: 
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with  x'  evaluated  at  the  stationary  point  x  on  Sj  and  z^  and  x^  evaluated  so 
as  to  make  the  phase  stationary. 


From  this  result,  we  see  that  the  determinant  is  positive,  so  that  the 


w'*  k"-  p'*  .*■**. 


eigenvalues  of  each  sign  must  occur  in  pairs.  Thus,  the  only  choices  for 
sig(i^jj)  are  -4  and  0  and  the  only  effect  that  the  signature  factor  can 
have  on  the  final  result  in  (31,32)  is  a  multiplication  by  -1  or  +1 , 
respectively.  In  Appendix  D,  J  show  that,  in  fact,  the  signature  is  zero 
and  the  multiplier  is  +1.  In  this  case,  the  integral  I(x)  defined  by  (32) 
be  come  s 


Kx) 


F(u))  exp  {ii>J^(x,x',x^,x^)}  dto  . 


(35) 


I  will  assume  that  the  original  source  was  impulsive.  Thus,  from  the 
assumptions  about  F((ij)  in  Section  2,  it  can  be  seen  that  I(x)  is  a  band 
limited  Dirac  delta  function  of  the  argument,  Hz ,  x ' ,  x^,  x^)  .  Therefore  I 
se  t 

I(?)  =  6„( 4 (x, X ' , X  , X  ) ) ,  (36) 

B  -  -  -s  -r 

where  I  have  used  the  subscript  B  to  remind  us  that  this  is  a  band  limited 
delta  function. 


The  function,  4,  is  equal  to  zero  on  the  surface  Sj .  Thus,  the  support 
of  this  delta  function  includes  ,  In  fact,  this  is  the  only  zero  and  it 
is  isolated.  To  see  why  this  is  so,  let  us  take  the  gradient  of  4  with 
respect  to  x,  with  x',x^  and  x^  defined  by  the  stationarity  conditions  (Cl): 
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In  this  equation,  each  sum  on  k  i s  zero  by  the  stationarity  conditions  (Cl). 


Thus>  the  total  derivative  with  respect  to  xj  is  just  the  partial  derivative 
with  respect  to  the  explicit  zj  in  #.  ffe  can  now  conclude  that  Vi  is  not 
zero  by  the  same  assumptions  as  were  made  to  conclude  that  V'i  (19)  was  not 
zero.  Consequently,  the  only  zero  of  i,  subject  to  the  stationarity 
conditions,  (Cl),  is  the  surface  Sj,  itself.  By  standard  rules  about  delta 
functions  we  can  now  write  Kx)  is  terms  of  a  delta  function  of  a  single 
arclength  variable  having  the  property  that  it  measures  arclength  along  a 
curve  normal  to  Sj .  If  we  denote  that  arclength  by  s,  then 


6  (s)  8  (s) 

lU)  =  — - - = -  .  (38) 

|Vi|  |VT(x.g  )  +  Vt(x.x  )| 


This  delta  function,  with  support  on  Sj  is  the  singular  function  of  the 
surface,  Sj ,  introduced  below  equation  (14).  Below,  we  will  denote  this 
function  by  Yj(x)*  Determination  of  the  singular  function  of  a  surface 
constitutes  mathematical  imaging  of  the  surface.  A  plot  of  the  band  limited 
delta  function  YjgCz)  will,  indeed,  depict  the  surface.  In  fact,  standard 
seismic  output  depicts  the  reflectors  by  plotting  their  singular  functions 
within  a  scale  factor.  By  using  the  result  (38)  in  (30)  with  6b(s)  replaced 
by  YjB(x)  one  obtains 


3a.  c  (x)  A(x',x  )  A(x',x  )  |  h(z.^)| 

- 1 - - - — - - 

3n  ^  c*(x')  A(x,x^)  A(x,x^)  | 


n*[V'x(x',x  )  +  V'x(x',x  )] 

_ ^  ^  r 

|V'x(x'.x^)  +  V'x(x',x^)|  * 


(39) 


Again  the  reader  is  reminded  that  x',  x,  and  x.  u^e  determined  here  as 


functions  of  x  by  the  stationarity  conditions  (Cl),  so  that  the  entire 


result  is  a  function  of  x.  At  this  point  we  have  a  result  which  images  the 


j  th  reflector  through  the  dependence  of  da^/dn  on  the  function 


Only  remains  to  determine  the  peak  amplitude  of  this  result  when  x  is  on  the 


reflector. 


Peak_^2l  itud^ 


I  will  now  analyze  the  multiplier  of  YBj(x)  tbe  peak  of  this 


function,  that  is,  on  the  reflector,  itself.  To  do  so,  let  us  first 


introduce  the  acute  angle  6  between  the  upward  normal  to  the  surface  and  the 


incident  and  reflected  rays  on  the  surface.  Note  that  the  downward 


gradients  V't(x',Xj)  and  V'x(x' ,x^)  make  angles  of  n  -  8  with  this  normal  and 


make  an  angle  of  20  with  one  another.  Therefore, 


n  •[V't(x',Xj)  +  V't(x',Xj.)] 


2co s  0 


c(x ') 


V't(x',x  )  +V't(x',x  )|*  =  - - -  1 1  +  cos  20  ] 

®  I  c  (z') 


p  CO  S  0p 

lc(x')  J 


Finally,  in  Appendix  E,  I  show  that 


,  ,  2  CO  s  0 

)  =  — rT~  •  *  on  s  . 

r  '  c (x)  ~  J 


By  inserting  the  results  (40,  41)  into  (38),  one  obtains  the  following 


result  for  9aj/0n  at  its  peak;  that  is  for  x  on  Sj : 
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(42) 


—  -  Aj  YjB(x)  .  X  on  S.. 

Vith  Aj  defined  by  (24),  this  result  confims  the  discussion  at  the  end  of 
Section  2.  Nanely,  the  modification  (16)  of  Beylkin' s  original  inversion 
equation  (4)  asymptotically  produces  a  singular  function  at  each  reflector 
multiplied  by  the  upward  jump  in  a(x)  across  the  reflector. 


I  will  now  proceed  to  evaluate  yjg(x).  itself,  on  Sj .  To  do  so,  I  use 
(35),  (38)  and  (40)  with  x'  =  x  to  conclude  that 


7E 


2cos  9 


1 

2n 


I  F(u])  d(i)  . 


(43) 


We  see  here  that  the  actual  numerical  value  at  the  peak  depends  on  the 

opening  angle  6  between  the  normal  and  each  of  the  rays  from  x  to  x  on  S- 

—  s  —  j 

and  x^  to  x  and  Sj.  Onf ortunately,  we  do  not  know  this  angle.  From  (41) 
and  (11),  we  recognize  the  first  fraction  in  (43)  as  just  the  factor 
relating  |u)|  and  |k|  in  our  heuristic  transformation  (eq.  11)  from  (u,^)  to 
k.  The  reason  for  the  appearance  of  this  factor  now  becomes  more  apparent. 
In  (43),  I  have  expressed  our  answer  in  terms  of  a  band  limited  frequency 
domain  delta  function.  I  have  done  this,  because  that  is  the  way  the  data 
is  provided  in  the  seismic  experiment.  On  the  other  hand,  rg(x)  should  be 
expressed  as  a  band  limited  wave  number  domain  delta  function.  This  first 
factor  is  the  scale  between  these  two  transform  variables. 

With  (40)  as  a  guide,  the  modification  of  our  processing  formula  to 


obtain  a  result  which  is  G-independent  at  the  peak  is  apparent.  f.e  need 


only  modify  the  integral  formula  in  (16)  by  eliminating  the  factor  of 
|Vt(x,Xj)  +  VtCx.Xj.)]!  This  effectively  changes  the  scale  of  the  normal 
derivative  so  that  we  are  no  longer  differentiating  with  respect  to  n,  but 
with  respect  to  i,(z)  =  #(x. x ' . x .. x, )  evaluated  at  the  stationary  point 

J  -  _  _  _5 

related  to  Sj ,  Thus,  symbolically,  I  set 

d  4  - 

A(x,x^)A(x,Xj.) 

(44) 

I  ioj  dw  F ' !i))  exp{-i(i)  IVtCx.x^)  +  Vt(x,x^)))  . 


da 

a? 


c*(x) 


8n 


The  symbolic  differentiation  here  is  with  respect  to  in  the 

neighborhood  of  Sj.  (No  such  differentiation  is  actually  performed  since  I 
only  mean  to  give  a  suggestive  name  to  the  indicated  operations  on  the  right 
side  of  the  equation.)  All  of  the  asymptotic  analysis  carried  out  above  for 
da/dn  is  identical  for  da/d4,  except  that  in  the  evaluation  at  the 
stationary  point  one  must  eliminate  one  factor  of 

|Vr(x,Xj)  +  Vt(x.x,)|  =  2cos9/c(x)  (eq.  (42)).  Consequently,  the  peak  value 
of  the  reflectivity  function  daj/dfj  in  response  to  data  from  the  j  th 
surface  is 


F(u))  diu 


X 


on 


(45) 


and  this  replaces  (43).  Since  we  know  the  filter,  we  know  the  integral  on 


the  right.  Correcting  the  peak  output  by  this  scale  factor  provides  an 
estimate  for  Aj.  In  (IS),  the  relationship  between  Aj  and  the  change  in 
sound  speed  across  the  reflector  is  stated.  Thus,  if  we  know  the  sound 
speed  above  the  reflector  we  obtain  an  estimate  of  the  sound  speed  below  the 
reflector  consistent  with  the  Born  approximation  on  which  the  analysis  was 
ba  sed. 

In  suBimary,  then,  I  have  proposed  a  modification  of  Beylkin's 
fundamental  inversion  formula  forto  Born-approximate  data  and  then  proceeded 
to  analyze  the  output  as  applied  to  such  data  by  asymptotic  methods.  My 
conclusions  are  as  follows. 

( i)  The  output  da/dn  is  proportional  to  the  sum  of  bandlimited  scaled 
singular  functions  of  the  reflectors  plus  possible  lower  order 
(smaller)  asymptotic  contributions. 

( ii)  At  the  peak  value  of  the  output  on  a  given  reflector,  the  scale  factor 
of  each  singular  function  is  the  peak  value  of  the  bandlimited  singular 
function  multiplied  by  the  jump  in  a  across  that  reflector.  The  peak 
value  of  the  singular  function  is  proportional  to  the  cosine  of  an 
unknown  angle.  By  a  slight  modification  of  the  inversion  operator,  the 
data  can  be  processed  for  a  function  I  defined  as  da/di.  whose  peak 
value  is  just  the  jump  in  a  multiplied  by  the  area  under  the  frequency 
domain  filter  of  the  data. 

I  repeat  that  F(<i))  is  typically  a  tapered  version  of  the  original 
source.  Thus,  this  integral,  within  smoothing  is  just  the  source  in  the 


time  domain,  F(t),  evaluated  at  t  =  0. 


I  return  now  to  the  question  of  a  discontinuous  background  sound  speed. 
Let  us  consider  the  effect  on  the  output  of  a  c(x)  which  might  be 
discontinuous  above  S.  but  remains  continuous  in  the  neighborhood  of  S;.  As 

«/  J 

long  as  the  ray- the  ore  ti  c  phases  and  amplitudes  used  in  our  inversion 
formulas  include  the  effects  of  those  disontinuities  —  e.g.,  refraction  of 
rays  and  transmission  coefficients  —  then  the  results  obtained  here  remain 
valid.  After  the  surface  Sj  has  been  identified,  the  effects  {f  that 
surface  are  incorporated  into  the  integral  operator  to  invert  for  points 
be  1 ow  S j  . 


m 


More  generally,  given  a  set  of  surfaces  of  discontinuity  for  the 
background  velocity,  the  upper  velocity  is  used  for  a  small  region  below 
each  input  surface  and  then  effects  of  that  surface  are  included  to  process 
deeper.  This  was  the  method  successfully  used  by  Docherty  [1985]. 
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5.  ASnPIOTlC  ANALYSIS  FOB  DB(aHf»V-APFBOXIIIATE  DATA 


Let  ns  considerj  now,  the  consequences  of  applying  the  inversion 
foimnlas  (16)  to  Kirchhof f-approximate  data  D(ai,§^)  f or  a  single  reflector, 
rather  than  to  Born-approximate  data.  On  the  one  hand,  such  an  application 
is  suspect,  since  the  derivation  of  the  results  (16)  was  based  on  a  Born 
approximation  of  the  solution  to  the  forward  scattering  problem  for  the 
upward  scattered  data  from  the  subsurface.  On  the  other  hand,  Kirchhoff 
data  is  not  constrained  to  small  perturbations,  but  only  to  high  frequency, 
which  we  have  used  throughout,  in  any  case.  Furthermore,  Kirchhoff  data 
more  accurately  represents  the  upward  scattered  field  from  a  single 
reflector  (especially  for  separated  source  and  receiver)  than  does  Born 
data.  Thus,  a  useful  output  of  this  analysis  will  allow  us  to  dispense  with 
the  small  perturbation  constraint  in  the  forward  scattering  problem,  a 
constraint  which  was  imposed  on  the  original  derivation  of  (16)  while 
providing  us  a  better  estimate  of  the  effect  of  applying  our  inversion 
operator  to  field  data. 

In  order  that  the  reflector  in  question  be  properly  located,  it  will 
still  be  necessary  that  the  background  sound  speed  not  vary  too  much  from 
the  "true"  sound  speed  above  the  reflector  in  question.  However,  this  is  a 
somewhat  less  severe  constraint  when  one  contemplates  a  theory  that  will 
produce  something  better  than  a  linear  approximation  to  the  velocity 
variations  and  we  could  then  contemplate  applying  a  three  dimensional  "layer 
stripping"  method  to  proceed  from  one  reflector  to  another,  progressively 


deeper  in  the  subsurface. 


Kirchhof f-approsimate  data  for  a  separated  source  and  receiver  and 
single  reflector  can  be  found  in  many  sources,  including  Bleistein  [1986], 


eq.  (49).  In  our  present  notation  the  result  is 


0(^,0))  ~  ioj  I  R(x',x  )  A(x',x  )  A(x',x  )  exp  {iw[T(x '  ,x.)  +  t(x',x^)]} 
I  o  o  I  s  r 

»  o 
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n-[V'T(x',x  )  +  V't(x',x  ))  dS'  , 

-  -s  -  -r 


In  this  equation,  R(x’,x,)  is  the  geometrical  optics  reflection  coefficient. 


R(x',x^)  = 


|3t(x' , x^) /3n  -  Jv  (x'+)  -  v  (x'-)  +  (x ' , x^) / 3nj 
jaT(x  ' ,  x^) /an  I  +  V  {x'+)  -  V  (x'-)  +  pt  ( X  ’ ,  x^) /anj 


,  (47) 


This  result  into  (16)  to  obtain  the  following  multifold  integral 
representation  of  the  output  da/dn: 


c*(x)  ^  lh(x.^)l  |Vx(x,i  )  +  Vt(x,x  )|  j 

- - —  d  5  ^ - - —  w  dll)  F(u)) 

8n  A(i,x  )A(i.x  ) 

1  J  Q  s  r 


R(i', 

*  c* 


X  )  A(x',x  )  A(i’.i  )  exp  (ill)  f(x,x',x  ,x  )} 

_s  -  -  _ 


n-[V'T(x',x  )  +  V’t(x',x  )]  dS'  . 

-  - s  -  -r 


This  result  is  to  be  compared  to  the  integral  (26)  for  the  analysis  of 
Born-approximate  data  from  a  single  reflector  S^.  We  can  see  that  the 


N'.-.V-V 
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at  the  peak.  By  using  this  observation  in  (43)  one  concludes  that 


da 

dn 


8R(x,x  )  cos*0  1 
-  -s _  _ 

clx)  2n 


F((i))  du) 


X  on  S 


(52) 


I  now  show  how  this  peak  value  is  related  to  the  previous  result.  (43). 
To  do  so.  I  use  the  fact  that  Vr  makes  an  angle  of  n  -  6  with  the  normal  at 
the  stationary  point  to  rewrite  (47)  as 


(53) 

-1  I  -a  -a  -a  a 

C  (x)  COS0  +  Jv  (x  +  )  -  V  (x-)  +  c  (x)  cos  0 

Now  let  us  assume  for  a  moment  that  the  jump  in  v(x)  across  S  is  small.  One 
then  obtains  the  leading  order  approximation  of  R(x,x^)  fox  small  values  of 


the  j  ump  by 

expanding 

the 

square 

root 

both  in  the 

numerator 

and 

the 

denominator , 

retai ning 

two 

terms 

in 

the  numerator 

and  one 

in 

the 

denominator.  The  result  is  that 


c  (x)  COS0  -  Jv  (x  +  )  -  V  (x-)  +  c  (x)  cos  0 


4R(x.x  )  cos  0  »  c  (x) 


v*(x-)  v*(x+)  J 


r(x+)  -  v(x~) 


2c  (x) 


=  A. 


(54) 


where  I  have  used  (15)  and  (24)  with  x  =  x'  to  obtain  the  1 


That  is.  to  leading  order 
the  previous  result. 


for  small 

f 


jumps,  the  present  result 
result  (52)  predicts  that 


ast  result, 
agrees  with 
the  output 


! 


H  ow  ev  e  r 


will,  in  fact,  produce  the  geometrical  optics  reflection  coefficient  at  the 
peak,  even  though  this  coefficient  is  a  nonlinear  function  of  the  jump  in 
sound  speed. 


On  the  other  hand,  an  estimate  of  the  angularly  dependent  reflection 
coefficient  is  not  useful  unless  we  have  a  means  of  estimating  the  angle. 
Fortunately,  we  do  have  such  a  means.  What  one  must  do  is  compute  daj/di 
[Eq.  (44)]  as  well  as  da/dn,  since  we  know  that  the  peak  output  of  the 
former  by  a  factor  of  2cos0/c(x).  That  is.  if  the  surface  S  is  one  of  the 
surfaces  Sj  then  the  peak  value  (45)  is  replaced  by  the  result. 


da.  ,  1 

— -  =  4R(x,x  )  cos  6  —  F(a))  dw  ,  x  on  S  . 

df.  "  ~  2n 

J 


Now,  by  taking  the  ratio  of  the  results  in  (52)  and  (55)  one  obtains 


da  . 

_ J 

dn 


2 cos  0 

— — t —  ,  X  On  S  =  S .  . 
c(x)  -  J 


(56) 


From  this  ratio,  cos&  is  computed.  Given  cos0.  (53)  is  used  to  compute  the 

sound  speed  below  the  reflector  in  terms  of  the  sound  speed  above  the 

reflector.  This  is  a  fully  nonlinear  estimate  of  sound  speed  consistent 

with  the  high  frequency  asymptotic  analysis  that  has  been  carried  out. 


For  esthetic  reasons  it  would  also  be  desirable  to  obtain  a  result  in 
which  the  peak  value  of  the  output  was  equal  to  the  reflection  coefficient 
multiplied  by  the  area  under  the  filter.  Py  examining  (52),  we  see  that  we 


need  only  eliminate  a  factor  of  8co8*0/c(x)  at  the  peak.  From  (50).  we  see 
that  the  way  to  do  this  is  to  introduce  into  the  inversion  operator  (16)  a 
divisor  of  c* (x) |Vt (x.x^)  +  Vr(x,x^)|*.  Thus,  I  propose  one  other  operator, 
which  I  call  the  reflectivity  function  and  denote  by  P(x) : 


i(ii  d(i)  F(<>))  exp {-i<ii[T(x,Xg)  +  ^(x,*^)))  D(£,(i))  . 


(57) 


With  no  further  analysis,  we  can  be  assured  that  the  peak  value  of  the 
reflectivity  function  is  given  by 

p(x)  ~  R(x,x^)  ^  I  F(w)  d(i),  X  on  S  =  .  (58) 

I  note,  further,  from  comparing  this  result  with  (55),  that  we  could  as 
easily  determine  cos0  through  the  ratio 

da . 

- ^  ~  cos  e,  X  on  S  =  S  .  (59) 

p(x)  ^ 

In  summary,  then,  the  operator  derived  on  the  basis  of  the  Born 
approximation  produces  a  reflector  map  when  applied  to  Kirchhoff  approximate 
data  with  a  nonlinear  estimate  of  the  jump  in  sound  speed  across  reflectors, 
consistent  with  the  geometrical  optics  reflection  coefficient  of  each 
reflector.  The  angular  dependence  of  the  reflection  coefficient  is  resolved 
by  simultaneously  computing  two  inversion  outputs  and  taking  their  ratio. 
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6 .  OONCLDSIONS 


Starting  from  an  inversion  operator  proposed  by  Beylkin,  I  have 
proposed  three  modifications  of  that  operator.  Each  of  those  operators  art- 
shown  by  asymptotic  analysis  to  produce  a  reflector  map  when  applied  either 
to  Born-approximate  input  data  or  Kirchhof f-approximate  input  data.  The 

peak  value  of  the  output  of  these  operators  is  proportional  to  the  Born 

approximate  reflection  coefficient  (or  jump  in  a(x))  in  the  former  case  and 
to  the  geometrical  optics  reflection  coefficient  in  the  latter  case.  The 

output  also  depends  upon  the  opening  angle  between  specular  rays  from  the 
source  to  the  reflecting  surface  and  from  the  receiver  to  the  reflecting 
surface.  This  opening  angle  is  unknown,  but  can  be  eliminated  for  either 
type  of  data.  These  results  are  valid  for  three  source/receiver 
configurations  of  interest:  common  (or  fixed)  source,  common  receiver  or 
common  offset,  with  the  last  of  these  including  the  zero  offset  or 

backscatter  case.  The  analysis  also  assumes  that  there  are  no  caustics  in 
the  ray  trajectories  from  subsurface  points  to  sources  or  receivers.  The 


;-i 


.% , 
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analysis  allows  for  a  curved  datum  surface. 
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FI6DBE  CAPTION 


Figure  Cl:  Triple  x',  x^,  x^  satisfying  stationary  phase  conditions, 

Eq.  (Cl)  for  a  horizontal  observation  surface,  horizontal  reflector  and 
constant  background  sound  speed. 
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APPENDIX  A:  KAT-IBEORETIC  GREEN’S  FONCTION 


In  this  appendix  I  present  some  results  about  the  WKBJ  or  ray-theoretic 
Green's  function, 

GCx.x'.taj)  =  A  (x ,  X  '  )e  xp  {iuiT  (x,  X  ' ) )  >  (Al) 

satisfying  the  inhomgeneous  Helmholtz  equation. 

2 

v’g  +  ^  G  =  6(x  -  x') .  {A2) 

c 

In  this  equation,  c  denotes  the  background  sound  speed,  c(x)>  with  its 
argument  omitted  when  it  is  clear  in  context.  The  travel  time  t  and  the 
amplitude  A  are  chosen  so  that  this  equation  is  satisfied  to  order  w*  and  la. 
The  equations  they  satisfy  are  called,  respectively,  the  eikonal  and 
transport  equations: 

IVxJ*  =  1/c*,  2Vt'TA  +  AV*x  =  0  .  {A3) 

A  solution  of  these  equations  is  required  to  describe  a  wave  emanating 
from  the  source  point,  x'.  Thus,  the  initial  conditions  for  t  and  A  are 

t(x'.x')  =0,  |x  -  x’|A(x,x’)  ->l/4n,  as  x  “>x'  •  {A4) 

The  problem  for  t  is  solved  by  the  method  of  charac ter i st icf •  [Bleistein, 
1984  ],  The  characteristics  satisfy  the  system  of  ordinary  differential 


equations. 


—  -2,  X  =  x‘ ,  for  0=0;  —  =  V[l/c  ] ,  £  =  Vx  . 

do  do 


Let  ns  denote  the  initial  valne  of  f  by  that  is> 


Po  =  pU'.s')  . 


This  initial  value  of  fi  oot  deteimined,  except  for  its  magnitude  which 
must  equal  l/c(x')  from  the  eikonal  equation.  Each  choice  of  the  direction 

of  p  or.  equivalently,  each  choice  of  p^#  and  p»o  *ithin  the  constraint  on 

|p|>  determines  or  '‘labels"  an  initial  direction  of  a  ray.  which  is  the 

solution  curve  x(o).  or.  more  completely  l(o.Pi»>Pao^  *  rays  "sweep  out" 

a  volume  of  space  as  we  vary  the  three  ray  parameters.  If  rays  for 
different  choices  of  (pio<Pse)  cross  one  another  —  i.e.,  if  there 

are  no  caustics  —  then  away  from  x  =  x'.  there  is  a  one-to-one 

correspondence  between  values  of  x  and  values  of  the  triple  (o.p^, .p,,).  It 
is  assumed  throughout  the  analysis  that  there  are  no  caustics  in  the  region 
of  interest. 


The  travel  time  also  satsifies  an  ordinary  differential  equation  with 
respect  to  o  —  i.e..  an  ordinary  differential  equation  along  the  ray  — 
namely. 


and  thus  is  determined  by  integration  along  the  ray.  The  transport 
equation,  the  second  equation  in  (A3),  can  also  be  written  as  an  ordinary 


*.  •*.  **1 
W  •  a 

•  >  «  a  ..  m  k 


• .  • .  • « » 
•  •  a  ■*  V  •  •  ■ 

^  ^  * 


A2 


differential  equation  along  the  ray>  with  solution. 


A(x,x')  = 


'♦"IPjoJ  I 


i/a  • 


The  Jacobi  determinant.  J.  is  given  by 


(AS) 


3(Xj^,  Xj.Xj) 

J  =  - - —  =  det 


The  assumption  that  there  are  no  caustics  assures  us  that  J  does  not  vanish, 
except  at  o  =  0,  that  is.  at  x  =  x'. 

There  are  two  practical  ways  to  calculate  J.  First.  J  can  be  related 
to  the  spreading  of  ray  tubes.  That  is,  J  is  proportional  to  the  normal 
cross-sectional  area  of  the  differential  ray  tube.  Thus,  if  one  is 
•shooting"  rays,  a  reasonably  accurate  approximation  of  J  can  be  determined 
by  measuring  the  area  of  the  tube  of  neighboring  rays  at  the  arrival  point. 
The  constant  of  proportionality  depends  upon  the  ray  parameters  that  are 
used  in  the  shooting  method. 

Alternatively,  the  components  of  the  Jacobi  matrix  from  which  J  is 
computed  also  satisfy  ordinary  differential  equations  with  respect  to  o; 
that  is  ordinary  differential  equations  on  the  rays.  See,  fcerveny,  etal. 
[1977]  or  Bleistein,  [1984], 


P 

dx 


^Pxo 

dx 

dp 

^  A 


(A9) 


One  can  see  from  this  brief  discussion  that  ray  theory  lends  itself 
naturally  to  initial  value  problems,  that  is,  to  the  'shooting  of  rays*  from 
a  prescribed  point.  On  the  other  hand,  in  the  application  to  inverse 

problems,  one  is  constantly  thinking  in  terns  of  a  travel  time  and  amplitude 
between  two  fired  points,  say,  r  and  x^,  for  example.  Indeed,  it  will,  at 
times,  be  convenient  to  think  of  the  rays  as  emanating  at  x  and  arriving  at 
and  sometimes,  the  reverse.  There  are  both  theoretical  and 

computational  difficulties  associated  with  the  transition  between  initial 
value  problems  from  each  point  and  the  boundary  value  problem  between  the 
two  points.  We  will  see  some  of  this  in  the  next  appendix. 


A4  - 


APPENDIX  B:  ANALYSIS  OF  h(x.^) 


The  purpose  of  this  appendix  is  to  introduce  simplifications  of  the 
determinant  defined  by  (6),  for  four  cases  of  interest  in  seismic 

exploration.  In  particular.  It  will  be  shown  how  this  matrix  is  related  to 
the  consitutents  of  the  ray  theoretic  Green's  function  developed  in  the 
previous  appendix. 

Casel  j _ Zero-Offset 

Let  us  suppose  that  we  are  considering  the  idealized  case  of  coincident 
source  and  receiver.  In  this  case. 


(Bl) 


the  sum  of  travel  times  becomes  just  twice  the  travel  time  to  the 
sour ce/re ce iver  point  and  the  determinant  in  (6)  becomes  the  simpler  result. 


h(x,^)  =  8  det 


Vt(x.x  ) 

-  "S 


Vr(x.x  ) 

~  “S 


Vt(x,x  ) 

ac  -  -s 


(B2) 


The  gradient  ■Vt(x,X5)  is  the  value  of  the  p-vector  (introduced  in  Appendix 
A)  on  the  ray  initiated  at  x^  and  arriving  at  depth  at  the  point,  x.  Since 
this  vector  satisfies  the  eikonal  equation  with  right  side  l/c*(x). 
independent  of  t,  we  have  the  result. 


To  exploit  this  identity  I  multiply  the  matrix  in  (B2)  by  another  matrix. 


(B4) 


before  computing  the  determinant.  Since  det  P  =  p, ,  the  result  of  this 
calculation  is 


h{x,^)  = 


8 

det 

Pt 

ap. 

P» 

ap. 

1 

n 

8  a(pj,Pj) 

PjC*(x) 

V 

p,c*(x)  0(5j,5g) 

ap. 

f\ 

u 

[uation. 

the  factor  1/c 

*(x) 

ari  se  s 

as  the  coefficient  in  < 

(B5) 


row  third  column,  which  makes  it  a  factor  of  the  entire  third  column. 


Let  ns  now  view  the  ray  connecting  x  and  Xj  as  starting  from  the  point 
at  depth,  x,  and  arriving  at  the  surface  at  the  point,  x^.  The  vector 
denoted  by  p,  here,  is  just  the  negatiy®  of  the  initial  value  of  the 
p-vector  on  that  ray,  say,  p*.  The  superscript,  s,  to  denote  a  ray  which 
emanates  from  x  and  arrives  at  Xg*  (Later,  1  will  have  need  of  the  notation 


•  \  ■-  V  '  ^rr  »r»  k-w  vna 


2^,  as  well.)  In  terms  of  this  new  vector,  the  result  (5)  can  be  rewritten 
as 


h(x,£)  =  - 


a(p  .p  ) 

*^10 


p,„c  (x) 


(B6) 


The  determinant  h(x,^)  has  been  written  in  terms  of  derivatives  of  the  two 
ray  parameters  at  depth  which  label  the  ray  propagating  from  x  to  Xg>  In 
fact,  under  the  assumption  that  there  are  no  caustics,  the  variables  ^  serve 
equally  well  as  parameters  to  label  the  ray  from  x  which  emerges  at  the 
upper  surface  at  the  point  *g(£) .  Thus.  h(x,^)  can  be  interpreted  as  being 
proportional  to  the  Jacobian  of  transformation  between  two  "ray-labeling" 
sets  of  parameters.  Each  parameter  pair  can  be  supplemented  with  a  third 
running  parameter  along  a  fixed  ray,  namely,  a,  of  Appendix  A.  Therefore, 
without  changing  the  value  of  h(x,£)  ,  two-by-two  determinant  in  (6)  can  be 
expanded  into  the  three-by- three  determinant. 


h(x.^) 


P..C 


8 

*(x) 


(B7) 


since  the  third  row  and  column  of  this  augmented  determinant  have  only 
zeroes  except  for  a  one  on  the  main  diagonal. 


Having  used  x  to  denote  the  initial  point  of  this  ray,  let  us  think  of 
the  running  coordinate  along  the  ray  as  being  some  other  variable,  say  1. 
The  derivatives  in  (B7)  are  to  be  evaluated  when  x  =  x_. 


Thus,  I  write 


3(Pio'P»o'®) 


J  |x  =  X3 


The  first  bracket  on  the  right  can  be  recognized  from  Appendix  A  as  jnst  the 
Jacobian  of  the  ray  family  initiated  from  x  and  evaluated  at  the  surface 
point,  Xg,  This  Jacobian  will  be  denoted  by  J(Xg,x)<  The  second 
determinant  can  be  re-interpreted  as  a  triple  scalar  product  of  two  tangents 
in  the  surface  S#  with  the  ray  tangent  dx^/da  =  p(x,,x).  The  cross  product 
of  the  tangents  is  in  the  direction  of  the  normal  vector  at  x^  and  has 
magnitude  ^g(x2)  =  ^8s'  complete  analogy  with  the  result  (29)  for  the 
surface  Sj .  Thus, 


a(Pi,.P,..«)  _  ®*s  ^  ®*s 

X  . 


lith  this  result,  h(x,£)  can  be  rewritten  as 


h(x,£)  =  - 


P^c*(x)  J(x^,x) 


(BIO) 


What  can  be  seen,  here,  is  that  the  matrix  h(x,£)  is  expressible  in  terms  of 
variables  which  are  computed  as  part  of  the  ray  theoretic  Green's  function 
and  in  terms  of  parameters  defining  the  structure  of  the  observation 
surface.  For  h  to  be  nonzero,  J  most  remain  finite  and  the  rays  should  not 
emerge  tangent  to  the  upper  surface.  For  h  to  remain  finite,  J  must  be 
nonzero  —  no  caustics.  The  singularity  at  p^,  =  0  is  only  apparent;  the 


product  PjjKXj/x)  must  remain  finite  in  the  limit  as  approaches  zero. 

This  product  is  exactly  what  appears  in  the  denominator  of  the  ray-theoretic 
amplitude.  (A8),  and  must  therefore  remain  nonzero  in  this  limit  because  one 
can  find  alternative  representations  of  the  amplitude  which  confirm  that  the 
amplitude  remains  finite  along  a  ray  that  is  initially  horizontal.  In  the 
special  case  where  the  upper  surface  is  flat,  the  dot  product  in  (BIO)  is 
simply  P,(Xg.x)  and  the  numerator  being  nonzero  requirf-s  that  the  ray 
direction  have  some  vertical  component  at  the  upper  surface.  Furthennore, 
in  this  case,  it  would  be  natural  to  take  f,  =  x,  and  I,  =  x,  and  then 
gg  =  1.  Consequently. 


•*) 


flat  upper  surface. 


(BID 


As  a  further  simplification,  if  c(x)  is  a  constant.  P|  is  a  constant  on 
each  ray.  In  this  case,  the  result  (Bll)  simplifies  still  further  to 


h(x,£)  =  -  ^c*J(x^,x)J  flat  upper  surface. 


constant  background. 


(B12) 


In  this  simplest  form,  the  reciprocal  dependence  between  h  and  J  provides 
the  total  structure  of  the  relationship  between  them. 


Case  2 :  Common  Source 


Let  us 


included  here  for  completeness. 


the  ca  se 

in  which 

the 

source 

point  is 

f  i  xe  d . 

Cohen. 

Blei stei n 

and 

Ragin 

[1985]  . 

It  is 

N'.'-'i-.v. 

For  this  case,  the  parametrization  of  the  source  and  receiver  points 


be  come  s 


=  const., 


(B13) 


with  the  function  ranging  over  the  observation  surface,  ,  as  J 

ranges  over  its  set  of  values  on  S^.  For  this  case,  the  determinant  h(z,{^) 


in  (d)  becomes 


l»(x,4)  = 


V[r(x,x  )  +  t(x,x  )J 
8  ““  r 


(B14) 


Thus,  while  the  first  row  remains  a  sum  of  two  vectors,  the  second  and  third 
row  are  identified  as  derivatives  of  single  vectors  exactly  of  the  form 
discussed  above.  Just  as  in  the  previous  case,  the  matrix  on  the  right  is 
to  be  multiplied  by  the  matrix  P  defined  in  (B4)  before  the  determinant  is 
c(»ipnted.  However,  I  now  have  to  be  more  precise  because  the  source  and 
receiver  are  separated.  The  p-vector  I  nse  in  P  is  the  vector  p(x  ,x).  By 


carrying  out  the  matrix  multiplication  and  then  taking  determinants,  the 
following  result  is  obtained: 

1  +  c*(x)p^.p^  I 


0 

In  this  equation. 


=  — 


det 


p  +  p 

is  Ir 


®Pir 


3^1 

3p 


ir 


p  +  p 


!!i£ 

^P»r 


I 

Ps  “  ®^-*-s^ ’  ®r  " 


(B16) 


Except  for  the  subscripts,  the  two-by-two  determinant  in  the  lower  left  hand 
corner  in  (615)  is  exactly  of  the  form  of  the  two-by-two  treated  in  Case  1. 
Ihus.  use  of  the  same  method  as  in  Case  1  yields 

1  +  c*(x)g*’  pj]  p(x  .x)*11  Jg 

- - - i - £ - £!-£  .  (617) 

J(Xr-S> 

The  subscript  zero  is  again  used  to  designate  vectors  at  x  pointing  along 
the  rays  to  x^  and  x^.  and  superscripts  on  the  p-vectors  to  denote  that  these 
are  the  vectors  on  the  rays  oriented  from  x  to  x^  and  Xj..  reverse  of  the 
orientation  on  the  subscripted  vectors  above.  The  vector  P(x^.x).  is  the 
final  value  of  the  vector  pj  at  the  surface  point.  Xj- 


h(x,^)  =  - 


As  above,  the  determinant  h(x.4)  is  expressed  in  terms  of  elementary 
parameters  of  the  ray- theoretic  Green's  function  and  parameters  which 


characterize  the  upper  surface.  Note  that  the  same  simplifications  as  were 
made  above  can  be  made  here,  too. 


When  this  result  is  compared  to  the  previous  result,  (BIO),  we  see  that 
there  is  one  additional  way  in  which  h  can  be  zero  in  (B17).  The  first 

factor  of  the  numerator  might  vanish.  This  can  only  occur  if  the  vectors  So 
and  pf  are  anti-colinear.  That  possibility  was  eliminated  in  Section  3. 

Case  3 :  Common Receiver 

For  the  common  receiver  case,  we  need  only  interchange  the  source  and 
receiver  point  in  the  previous  case.  When  this  is  done  in  (16),  the 
following  result  is  obtained: 

+  c*(x)p^*  p®]  p(x  ,x)*ajg~ 

- 1 - U - ! - .  (B18) 

p®,c*(x)  J(x^,x) 

Case  4 :  The  General  Case 

When  both  and  x^  are  more  general  functions  of  t,  there  is  not  a 
great  deal  of  simplification  that  can  be  achieved  in  the  representation  of 
h(x,^) .  The  determinant  h  in  (6)  can  be  rewritten  as  a  sum  of  four 
de  terminants  : 


- 


h(x,t)  =  h  (x,S)  +  h  (x.g,)  +  h,(x,£)  +  h Ax.i)  > 


(B19) 


VIt(x,x  )  +  x(x,x  )1 
-  -s  -  ~T 

It7 

I{7 


vtx(x,x  )  +  t(x,x  )J 

-  -s  r 

k 

VtT(x,X  )  +  t(x,x  )] 

-  -s  ~  r 

k 

V[x(x,x^)  +  x(x,x  )] 

-  r 


-lx-  Vx(x,x  ) 


entire  determinant  in  the  case  of  a  common  source;  hi(z,^)  is  the 

same  determinant  as  for  the  common  receiver  case.  Thns.  (B18)  and  (B17)  can 
be  used  to  write  these  determinants  in  terms  of  the  parameters  of  ray 
propagation: 


=  - 


(B24) 


=  - 


[l  +  c*(x)e*«  2^j  £(5^.*)-®,,  «r 


(B2S) 


The  remaining  two  determinants,  h,(X,£)  and  h^(X,£)»  are  not  so  easily 
dealt  with.  The  reason  is  that  the  second  and  third  rows  do  not  involve 
derivatives  of  the  same  travel  time.  That  is,  one  row  has  t(x,Xj)  while  the 
other  has  t(x,x^).  Thus,  a  simplification  to  merely  a  ray  determinant,  J, 
is  not  possible  since  these  rows  are  related  to  different  ray  families.  On 
the  other  hand,  each  of  the  matrix  entries  is  separately  expressible  in 
terms  of  ray  parameters  and  surface  parameters.  That  derivation  now 
follows. 


Let  us  consider  a  typical  term,  3pj(x,x)/a?j,  for  i, j  =  1,2,  and  x  to 
be  evaluated  along  the  ray  to  (and  ultimately  at)  or  z^..  It  is  not 
necessary  consider  derivatives  of  p,,  because  the  eikonal  equation  (A3)  can 
be  used  to  express  them  in  terms  of  p,  and  Pj* 


(B26) 


9p,  _  Pj  ap,  p,  ap, 

»?7  *  ■  ■  57  HI-  • 

For  the  derivatives  of  interest. 


ai,  a?. 


*  i»j  “  lj2, 


(B27) 


When  W  is  evaluated  at  or  xr,  the  derivatives  dx^/d^j  simply  become 

derivatives  of  the  expressions  Xj(£)  and  x,.(^)  with  respect  to  the 
components  of  5.  Thns,  these  derivatives  are  known  from  the  surface 

geometry.  For  the  other  factor  on  the  right  in  (B27),  I  again  take  the 

point  of  view  that  f(x,x)  is  the  negative  of  the  initial  value  of  the 

p-vector  along  the  ray  starting  from  x  and  propagating  to  denoted  above 
By  pI  or  pg.  As  noted  in  Appendix  A,  the  first  two  components  of  these 
vectors  may  be  used  as  ray^labeling  parameters:  see  the  discussion  below 
(A6)  .  Thus,  as  a  first  step. 


3P.  3p., 

1  1* 


(B28) 


Now,  as  in  Appendix  A,  the  family  of  rays  emanating  from  x  and  covering  a 
volume  around  that  point  are  viewed  as  defining  a  transformation  or  change 
of  variables  from  i  to  (pio'Pao *  with  the  inverse  of  that  transformation 
mapping  (Pig,Pjg,o)  to  x.  The  expression  in  (B28)  is  an  element  of  the 
Jacobi  matrix  of  the  latter  transformation.  On  the  other  hand,  the  matrix 
of  the  former  transformation  is  more  familiar.  It  is  the  matrix  whose 


Bll 


dcterainant  J(5,2)  plays  a  crocial  role  in  tlie  propagation  of  anplitnde 
along  rays.  Thns,  let  ns  introduce  the  aiatriz 


PijJ  ■  If..  .  j.t  [ly]  .  I(f.s) 


(B29) 


.  jm  . 
•••  N/;-’ 

• '  f.,- / 


Since  the  two  transformations  in  question  here  are  inverses  of  one  another, 
the  matrices  of  transformation  mast  also  be  inverses  of  one  another.  (In 
making  this  claim.  I  am  assuming  that  there  are  no  caustics  in  the  ray 
family  over  the  domain  of  interest.)  Urns,  the  derivative  in  (B28)  can  be 
rewritten  in  terms  of  the  inverse  of  the  matrix  in  (B29).  The  result  is 


(B30) 


!ii  = 

®Pi«  -1  ~ 

—  -  -  T  -r)  /.rtf 

a^k 

J  ^2*  =  '  vOI 

axk 

.  ®Pio 

lii' 

.V.'  •-  o  ■ 

»  k  <■  aS  .V  J 


In  this  equation,  cof  denotes  cofactor  in  the  matrix  tly]  in  (B29) 


This  result  is  now  used  in  (B26)  to  write 


i  -1  -  r  k  k 

=  J  (x,x)  2  9^  cof  -  '  i»j  = 


(B31) 


^  *  V  * 


B12 


4^ _ Two-and-one-half  dimensions 


There  is  a  whole  class  of  special  cases  where  can  be  written 

more  simply.  These  are  the  so~called  two  and  one  half  dimensional  (2.SD) 
cases.  Here,  it  is  assumed  that  the  earth  parameters,  c  and  a  are  functions 
of  (z,z)  only  (essentially  two  dimensional)  while  the  wave  propagation  is 
three  dimensional.  Furthermore,  only  one  line  of  data  is  gathered,  say,  for 
y^  =  y^  =  0.  However,  for  such  a  medium,  all  lines  parallel  to  the  line 
y  -  0  would,  of  necessity,  produce  identical  data  and  only  the  output  of  the 
algorithn  in-plane,  that  is,  for  y  =  0,  is  required. 


In  the  inversion  foimnla.  (Id),  then,  the  data,  D(£,<»),  is  independent 
of  and  the  integration  in  can  be  carried  out  by  the  method  of 
stationary  phase.  The  stationary  point  turns  out  to  be  ^  and  the 

remainder  of  the  integrand  need  only  be  evaluated  in-plane.  In  this  case, 
h(x,^)  can  be  evaluated  in  terms  of  the  in-plane  ray  Jacobians  associated 
with  two  dimensional  wave  propagation.  Out-of-plane  effects  are  accounted 
for  through  a  scaling  by  )|o.  In  particular,  h(z,£)  is  expressible  in  terms 
of  2X2  determinants  with  the  first  row  always  being  related  to  a  p-vector. 
Thus,  a  case  like  h,(x,£)  or  h4(x,^)  cannot  occur  and  h  is  expressible  in 

terms  of  the  two  Jacobians  associated  with  rays  between  x  and  x  or  between 

-  -s 

X  and  x^> 


Ji^ 

> 

^  “J 


.v; 


h.  ’-j  'jt. 


-*• 


In  Bleistein,  Cohen  and  Hagin,  [1985],  the  details  of  this  computation 
were  carried  out  for  the  cases  of  common  source,  common  receiver  and  common 
offset.  It  should  be  noted  that  the  2.SD  common  source  or  common  receiver 
cases  do  not  arise  directly  from  the  3D  common  source  or  common  receiver 
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case.  In  the  3D  common  source  case,  for  example,  there  is  only  one  source 
in  the  plane  with  y  =  0.  Thus,  the  data  gathered  on  lines  parellel  to  y  =  0 
is  not  identical.  To  make  lines  of  identical  data  as  the  2.SD  model 
requires,  the  source  must  be  moved  to  each  new  line,  y  =  const,  when  the 
receivers  are  moved.  Thus,  2.SD  common  source  corresponds  to  the  3D  case  in 
which  there  is  a  line  of  sources,  say  with  x  =  0,  and  the  data  from  each 
source  is  gathered  along  an  orthogonal  line,  y  =  const.  Now  each  line  of 
experiments  will  be  identical. 

The  results  for  h(z,^)  2.5D  will  not  be  restated  here  because  the 

entire  processing  formula  in  2.SD  changes.  That  result  is  an  integral  over 
the  source/receiver  line  (or  curve)  with  an  adjustment  of  the  inversion 
integrands  provided  in  this  paper  to  account  for  the  out-of-plane  stationary 


phase  computation. 


APFBNOIZ  C:  ANALYSIS  OF  TBE  STATIONAKI  MASE  CONDITIONS 


The  purpose  of  this  Appendix  is  to  discnss  the  conditions  of  stationary 
phase,  that  is,  the  conditions  that  the  fonr  first  partial  derivatives  of  i 
in  (30)  are  all  eqnal  to  zero.  The  reader  is  reminded  that  the  travel  time  t 
is  symmetric  in  its  initial  and  final  coordinates,  Thns,  each  of  the 
gradients  appearing  in  (30)  is  a  p-vector  directed  tangent  to  the  ray.  For 
example,  V^x{x',x^)  =  p(Xj,x'),  is  a  p-vector  tangent  to  the  ray  from  x'  to 
x^  (from  second  argument  to  first  argument)  evaluated  at  x^  (evaluated  at 
first  argument).  It  has  magnitude  I/c*(x^)  and  is  directed  away  from  the 
initial  point,  x'.  Similarly.  ^''c^i'.Xj)  =  p(x’»Xj)  is  evaluated  at  x',  has 
magnitude  l/c*(x')  and  is  directed  away  from  Xg> 


The  result  (30)  and  the  notation  for  gradients  introduced  here  are  used 
to  write  the  conditions  that  the  phase  be  stationary  as  follows: 


dx 


‘‘is 


'  d^  ^  *  w  ■  *  IT  '  dT  ' 


[  b(x-.x^)  +  p(x'.x^)  ]  •  If-  =  0  .  b-1,2. 


(Cl) 


It  is  assumed  that  a  proper  parameterization  has  been  used  for  which  the  two 
vectors  in  each  case  (m  =  1,2)  are  linearly  independent. 


The  second  condition  is  easier  to  interpret.  It  states  that  the 
tangents  to  the  rays  from  x^  and  x^,  to  the  surface  point  x’  have  equal 
projections  on  two  linearly  independent  tangents  in  the  reflecting  surface. 
Consequently,  the  projections  of  these  two  vectors  onto  Si  must  be  equal. 


This  is  jast  Snell's  law  for  reflection.  The  magnitndes  of  the  p-^rectors 
must  be  equal  (to  l/c*(x'))  and  hence  the  out-of-plane  components  must  be 
equal  in  magnitude,  as  well.  Indeed,  the  normal  components  of  these  vectors 
are  of  the  same  sign  and  must,  in  fact,  be  equal. 

The  first  condition  in  (Cl)  ties  the  points  on  the  two  surfaces  to  the 
output  point.  X.  Let  ns  consider  the  rays  from  x  to  the  upper  surface 
points,  x^  and  x^.  Similarly,  we  consider  the  rays  from  x'  to  the  upper 
surface  points,  x^  and  x,..  For  each  pair  of  rays  we  take  projections  on 
tangents  at  their  respective  emergence  points.  The  sum  of  these  projections 
for  each  pair  of  rays  must  be  equal  to  one  another.  This  must  be  true  for 
two  linearly  independent  tnagents  at  each  point. 

At  first  glance,  it  may  not  seem  apparent  that  such  a  condition  can 
over  be  satisfied.  However,  consider  the  case  in  which  i  is  on  the 

reflecting  surface.  S j ,  Then,  for  x'  =  x  (and  a  chosen  accordingly)  the  two 
pairs  of  rays  overlay  one  another  and  these  stationarity  conditions  are 
automatically  satisfied  for  any  pair  of  surface  points,  x^  and  x^..  Thus,  we 
would  only  have  to  find  such  a  pair  for  which  Snell’s  law  is  satisfied,  as 
well.  Indeed,  if  there  were  no  such  pair  in  the  seismic  experiment  being 
modeled,  then  that  subsurface  point  would  not  be  one  for  which  the 

stationarity  conditions  are  satisfied  a^  that  point  would  not  be  imaged. 

On  the  other  hand,  there  are  many  candidates  for  source /receiver  pairs 
on  the  upper  surface  when  x'  =  x*  To  find  them,  proceed  as  follows.  At 

x'  =  X,  pass  a  plane  through  the  normal  to  Sj.  In  the  plane,  choose  two 

directions  making  equal  angles  with  the  normal.  Dse  these  as  initial 


directions  for  rays  from  the  point.  Snell's  law  is  satisfied  for  this  pair 
of  rays.  The  pair  of  emergence  points  at  the  upper  surface  are  candidates 
for  a  source/receiver  pair.  Vary  the  opening  angle  of  the  ray  pair  in  the 
normal  plane  and  rotate  the  plane.  Thereby,  obtain  a  two~-diinensional 
continuum  of  candidate  source/receiver  pairs  in  the  upper  surface. 

Let  us  suppose  now  that  such  a  pair  is  available  in  a  given  seismic 
survey  when  x  is  on  Sj .  Given  that  pair,  it  is  argued  by  continuity  that 
for  X  nearby  Sj  there  must  by  points  x',  x^  and  x^  satisfying  (Cl)  and 

nearby  the  solution  obtained  in  the  limit  when  x  is  on  Sj. 

Constant  Background  Soundspeed 

Further  insight  into  the  stationarity  conditions  is  gained  by 
considering  the  case  of  constant  background  speed  and  flat  layers,  as  in 
Figure  Cl.  Given  a  point,  x,  a  perpendicular  is  dropped  to  the  surface  S  j . 
This  determines  a  point,  x'.  Pass  a  plane  through  the  normal  and  draw  the 
rays  at  equal  angles  to  the  upper  surface.  This  determines  a  pair  of 

points,  as  candidates  for  x^  and  Xj.*  Por  this  pair  of  points,  the  sum  of 

projections  on  either  side  of  the  first  line  of  (Cl)  is  equal  to  zero. 

Thus,  this  triple  of  points  satisfies  both  conditions  of  stationarity. 

The  three  points,  x',  x^  and  Xj.  must  be  in  the  same  plane  in  order  that 
Snell's  law  be  satisfied.  If  x  were  not  in  the  same  plane,  then  the 
projections  of  its  p-vectors  would  no  longer  be  colinear  and  could  not  siun 
to  zero.  On  the  other  hand,  the  sum  of  projections  of  the  p~vectors  from  x' 
would  remain  zero.  Thus,  the  first  conditio.-,  in  (Cl)  could  not  be 


satisfied.  Similarly,  if  x  is  in  the  normal  plane  but  not  on  the  normal 
line  the  first  condition  could  not  be  satisfied.  That  is,  the  conditions  of 
stationarity  are  satisfied  by  three  points  x*,  x^  and  x^.  which,  along  with  x 
lie  in  a  plane  normal  to  the  reflector  with  x*  at  the  foot  of  the  normal  to 
Sj  drawn  from  x>  only  freedom  left  in  these  conditions,  then,  are  the 

opening  angle  of  the  rays  at  x'  and  the  orientation  of  the  normal  plane. 
Below,  I  discuss  how  these  are  further  constrained  for  particular 
source /receiver  configurations  and  this  flat  reflector  constant  background 
model. 

Case  1:  Zero-Offset 

When  the  source /receiver  pair  are  coincident,  the  opening  angle  of  the 
rays  at  x*  must  both  be  zero;  both  rays  from  x'  to  x^  and  x^  must  be  the 
normal  ray  to  the  surface,  passing  through  x*  The  stationary  point  on  the 
upper  surface  and  the  point  x'  must  have  the  same  transverse  coordinates  as 
X,  itself.  The  stationarity  conditions  are  completely  satisfied  by  these 
three  points.  Because  of  the  degeneracy  of  this  case,  a  specific  normal 
plane  is  not  determined.  However,  that  is  secondary  to  determining  the 
actual  triple  of  points,  itself. 

The  generalization  of  this  result  to  curved  surfaces  and  variable 
background  is  fairly  straightforward.  Given  x.  find  a  normal  ray  from  Sj 
which  passes  through  x.  The  initial  point  of  that  ray  on  Sj  is  the  point 
x*.  The  point  where  the  ray  emerges  on  the  upper  surface  S,  ia  the 
source /receiver  point  which  completes  the  triple  of  points  satisfying  (Cl). 
For  X  on  Si>  there  is  clearly  only  one  stationary  triple.  On  the  other  hand 


for  X  on  the  evolute  of  S-  (the  envelope  of  normals  to  S^)  there  will  be 
more  than  one  triple.  In  order  for  the  asymptotic  methods  used  here  to  be 
valid,  it  is  necessary  to  assume  that  this  evolute  is  a  few  wavelengths  (at 
least  three)  away  from  Sj.  Thus,  it  is  assumed  that  the  reflector  is  not 
severely  curved;  that  is,  the  Btincipal  radii  of  cnrvatur£  of  the  reflector 
must  be  a  few  wavelengths  long. 

Case  2 j _ Common_ Sour c e 

Let  ns  suppose  now  that  the  source  point  is  fixed.  Given  x^  drop  the 
normal  to  Sj  and  thereby  determine  x'.  Pass  a  plane  through  x,,  x  and  x'* 
This  plane  is  normal  to  S.,  Draw  the  ray  from  x  to  x'.  Draw  the  reflected 
ray  in  the  given  normal  plane.  The  emergence  point  on  S,  is  the  point  X^* 
If  X  is  on  Sj,  set  x'  =  x  and  use  the  normal  at  that  point  and  the  fixed 
point  x^  to  determine  the  normal  plane.  Then  proceed  to  determine  x^  as 
before,  with  x  not  on  S j . 

In  a  theoretical  model,  receivers  are  spread  over  the  entire  upper 
surface.  In  practice,  the  spread  is  finite.  Thus,  the  spread  need  not 
extend  to  the  determined  x^.  In  that  case,  the  determined  point,  x'  will 
not  be  part  of  a  triple  satisfying  (Cl)  and  will  not  be  imaged.  In  the  text 
I  have  proceeded  as  if  such  candidate  points  are  indeed  stationary  points. 

As  above,  I  argue  by  continuity  that  for  curved  surfaces  and  variable 
c(x),  differing  'not  greatly'  from  the  constant  background  case,  the 
essential  features  of  this  analysis  still  apply. 


C5 
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Cate  rnwmon  Receive^ 

As  in  the  previous  appendix,  one  need  only  interchange  the  subscripts  s 
and  r  in  the  discussion  of  Case  2  tc  obtain  a  completely  analogous 
conclusion  here. 

S* ^ •  Common  Offset 

It  is  assumed  that  all  of  the  offset  pairs  lie  on  lines  that  are 
parallel.  We  rotate  the  normal  plane  containing  x  and  x'  until  it  is 
parallel  to  this  set  of  lines.  Indeed,  the  intersection  of  the  normal  plane 
and  the  upper  surface  contains  one  of  those  lines.  ."lOosc  the  opening  angle 
of  the  rays  from  x'  so  that  the  rays  emerge  ar  tiio  upper  surface  at  a 
separation  distance  equal  to  the  comiuon  offset  distance.  The  emergence 
points  are  the  pair  and  x^* 

S  Common  Midpoint 

There  will  only  be  a  solution  to  (Cl)  in  this  case  if  the  common 
midpoint  and  z  lie  along  a  coaunon  normal  to  Sj.  Furthermore,  in  that  case, 
all  source/receiver  pairs  are  stationary  points.  The  method  of  stationary 
phase  breaks  down  since  the  stationary  points  are  no  longer  isolated.  This 
is  a  case  which  requires  further  investigation. 
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APPENDIX  D:  NATBIX  SIGNATDBE 


The  purpose  of  this  appendix  is  to  show  that  the  signature  of  the 
matrix  defined  by  (33)  is  equal  to  zero.  To  do  so>  I  consider  first 

the  special  case  in  which  the  background  sound  speed  c  in  the  region  between 
the  upper  surface  and  the  reflecting  surface  is  constant,  the  layers  are 
flat  and  there  is  zero  offset  between  sources  and  receivers.  In  this  case, 
the  upper  surface  and  the  reflecting  surface  are  defined,  respectively,  by 


X 

Is 


X 


X  =  0, 
*r 


=  a 


x'  = 
a 


=  e 


(Dl) 


Furthermore,  the  travel  times  are  just  the  distances  between  initial  and 
final  point,  divided  by  c: 


=  Is  -  =  U  ~ 

t(5',Xjj)  =  (x*  -  xj/c.  t(x',x^)  =  lx'  -  x^|/c. 


(D2) 


These  results  are  used  to  simplify  as  defined  by  (9)  and  then  to 
computed  the  determinant  in  (33).  The  analysis  is  further  simplified  by 
setting  x'  =  x.  The  result  is 


0 

0 

-1/Hc 

0 


0  -l/Hc  0 

0  0  -1/Rc 

0  1/Hc  0 

-1/Hc  0  1/Hc  . 


(D3) 


For  this  matrix  it  is  fairly  straightforward  to  calculate  the  characteristic 
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equation.  The  result  is 


det[i^^  -  ^  ij  =  IX(1  -  X)  +  l]*/lHcJ*  =  0  ,  (D4) 

This  equation  has  two  double  roots.  X  =  [1  ^  \|5]/2.  Since  two  of  the  roots 
are  positive  and  two  are  negative,  =  0. 

Let  ns  now  consider  deforming  this  constant  background,  zero  offset, 
flat  layer  model  into  the  true  model.  If  the  signature  is  to  change  as  the 
model  is  deformed,  then  at  some  point  in  the  deformation,  at  least  one 
eigenvalue  must  be  zero.  In  fact,  exactly  two  eigenvalues  would  have  to  be 
zero  at  this  point,  since  det[$^g]  is  nonnegative  and  by  assumption,  the 
signature  changes. 

In  the  next  section,  it  is  shown  that  det[$^^]  is  proportional  to 
h(z,^).  It  has  been  assumed  that  h  is  nonzero  for  the  true  model.  1  now 
add  to  that  the  assumption  that  our  true  model  is  not  so  severely  different 
from  the  flat  earth  case  for  h  to  have  passed  through  a  zero  on  the  way  from 
one  model  to  the  other.  Thus,  sig[$^^]  =  0  for  the  true  model,  as  well. 
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A?PBND1X  B:  BELATION  BETVEBN  AND  detli^gl  AT  IME  STATIONARY  POINT 


In  this  appendix  (41)  will  be  verified.  To  do  so.  isnecessary  to 
evaluate  |h(x,4)  |  as  defined  by  (6)  subject  to  the  stationarity  conditions. 
(Cl)  and  the  additional  condition  that  x  =  x'.  As  a  first  step,  x  is 
replaced  by  x'  in  (6)  and  that  equation  is  rewritten  in  terms  of  p-vectors. 
The  result  is 


pCx'.x  )  +  p(x',x  ) 

-  -s  r 


Ji££ 

•  *  k  •  w 


ir'.V.V. 

f- 

■s." 

.  I 


!>{*', i)  =  det  fp(s'**j)  + 


To  calculate  this  determinant,  the  matrix  is  multiplied  by  a  matrix  whose 
determinant  is  known.  That  matrix  is 


dx'  dx'  dx 
do7  '  '  la 


where  each  vector  here  represents  a  column  of  K.  We  remark  that 


dx'  dx' 

det  K  =  n  •  —  X  —  =  Jg . 

dOj  do^ 


with  the  second  equality  being  equivalent  to  (29). 


Now,  in  multiplying  K  by  the  matrix  in  (El),  we  see  that  the  first  two 
elements  of  the  first  row  are  both  zero  by  (Cl),  while  the  third  element  is 
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^  m  r_i  : 


V'.'V  V -J 


given  by 


[efx'.Sg)  + 


n 


2cos  ft 
c(x) 


(E4) 


which  follows  from  (40).  Thus,  in  expanding  the  determinant  of  the  product 
by  the  first  row,  it  is  only  necessary  to  consider  the  lower  left  2X2  matrix 
after  multiplication.  Thus,  let  us  now  consider  a  typical  term. 


dz ' 

m  k  m 


a^j^ao 


m 


k,  m 


1,2. 


(E5) 


It  now  follows  that  if  the  matrix  in  (El)  is  multiplied  by  the  matrix  E 
before  calculating  the  determinant,  the  following  result  is  obtained; 


dot  h(x',^)  ^gj  = 


2 cos  0 
c  (x ' ) 


det 


(E6) 


for  x'  =  X  on  Sj.  The 
right  equality  in  (41) 


outer  equality  in  (41)  follows  from  this  result, 
follows  from  (40). 


The 
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ABSTRACT 


This  paper  is  based  on  a  paper  by  Beylkin  in  which  a  leading  order 
asymptotic  theory  for  inversion  of  acoustic  data  is  presented.  The  method 
is  based  on  earlier  work  by  Beylkin  in  the  theory  of  pseudo-differential 
operators  and  generalized  back  projections  or  Radon  transforms.  The  back 
projection  or  inversion  is  carried  out  with  respect  to  a  general  [c(x,y,z)J 
background  sound  speed.  The  asymptotic  limit  of  interest  is  high  frequency. 
The  inversion  operator  is  given  as  an  integral  of  the  observed  data  over 
frequency  and  over  the  observation  surface.  Beylkin  claims  that  his  result 
is  useful  for  finding  discontinuities  in  the  sound  speed,  but  be  does  not 
make  clear  how  this  is  to  be  done  in  practice.  I  show  how  to  modify 
Beylkin's  inversion  operator  to  obtain  an  operator  whose  output  is  an  array 
of  singular  functions,  one  for  each  reflector  (discontinuity  surface  of  the 
sound  speed)  in  the  subsurface.  The  singular  function  of  a  surface  is  a 
Dirac  delta  function  whose  support  lies  on  that  surface.  Thus,  the  array  of 
singular  functions  produces  a  reflector  map  of  the  subsurface.  The 
validity  of  modification  of  Beylkin's  inversion  operator  is  verified  by 
applying  it  to  band  limited  Born-approximate  and  then  Kirchhoff-approximate 
representations  of  the  upward  propagating  wave  field.  Multi-dimensional 
stationary  phase  is  applied  to  the  spatial  integration  over  the  variables  of 
the  field  representation  and  the  variables  of  the  observation  surface.  It 
is  confirmed  that  the  output  is  proportional  to  the  band  limited  singular 
functions  of  the  reflectors  and  further  that  one  can  estimate  the  jump  in 
velocity  across  each  reflector  from  the  peak  amplitude  of  the  output  on  each 
reflector.  This  is  done  for  the  cases  of  common  (or  single  fixed)  source, 
common  receiver,  and  common  (or  fixed)  offset  between  source  and  receiver, 
with  zero  offset  or  backscatter  as  a  special  case  of  the  last  of  these. 


